{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Predictive Control\n",
    "The general idea of figuring out what moves to make using optimisation at each time step has become very popular due to the fact that a general version can be programmed and made very user friendly so that the intricacies of multivariable control can be handled by a single program.\n",
    "\n",
    "In this notebook I will show how a single time step's move trajectory is calculated. We'll use the same system as we used for the [Dahlin controller](../6_Discrete_control_and_analysis/Dahlin%20controller.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import scipy.signal\n",
    "import scipy.optimize\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with a linear model of the system\n",
    "\n",
    "$$G=\\frac{1}{15s^2 + 8s + 1}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = scipy.signal.lti([1], [15, 8, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c146a8390>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPW9//HXJ5OdbEBChJAYkF1l0YgLVrAuxRXbq71ovXWra221tb3a3l7XX6+t2ltrpW63rq1Fa9VSi+ICuCsEQWQLhLAlbIEkELIn8/39MQONGMgAk5yZyfv5eMwjM2e+ybw9kLdfzpw5X3POISIisSXO6wAiIhJ+KncRkRikchcRiUEqdxGRGKRyFxGJQSp3EZEYpHIXEYlBKncRkRikchcRiUHxXr1wdna2Kyws9OrlRUSi0oIFC7Y553I6G+dZuRcWFlJcXOzVy4uIRCUzWxfKOB2WERGJQSp3EZEYpHIXEYlBKncRkRjUabmb2ZNmttXMluzjeTOzh8ys1MwWm9kx4Y8pIiIHIpSZ+9PA5P08fxYwNHi7Bnjk0GOJiMih6LTcnXPvAVX7GTIFeNYFfAJkmVn/cAUUEZEDF47z3POADe0elwe3bQrDzxYRCRvnHK1+R0ubn5Y2R2ubn1a/o7nVT5vf0eoPbA/cd7T5/bS2Be63f+x3ux//a6x/91f3r+2B++zZ5neBcaeNzGVMflaX/reGo9ytg20dLsxqZtcQOHRDQUFBGF5aRKJNa5uf+pY2GpqDt5bArbHd/YbmNhpb/TS1tNHY0kZTqz9wa3e/udVPU+uXHze3+mlu8wfKu9VPc9vuIvfvKfRI0C8jOSrKvRzIb/d4ILCxo4HOuceBxwGKiooiYy+LSEgaW9rY2djCzoYWdjS0UtvYws7GwNddja3UNrayqyl4a2ylrjlwv66plbqmNuqbW6lrbqO51X/Ar20GyfE+khLiSIqPIyneR2J84H5ifByJvjjSk+NJio8jwRfYluAL3JLi44iPM+J3bw/eT/AZCb444n1GQlzga7wvODYu8FxcnJEQZ/jijHif4YsLPB9nRoLPiIszfNbuefvXtrjg9wXus2ecWUfz4fALR7nPAG40s+nA8cAO55wOyYhEsMaWNqrqmtm+q5ltu5qoqmumur6Z7XXN1NQ3Bx+3sKO+hZqGZmrqW2jqpJTjDNKS4klPTqBXko+0pHjSkuLJTU8mNclHr8R4UpN8pCbEk5roIyXRR2rwlpzgIyUhsC05wUdyvI/kxDiSE3yBAvfFdVspxopOy93M/gJMArLNrBy4A0gAcM49CswEzgZKgXrgiq4KKyL719jSxqYdjWze0ciWnbtvTWytbWRrbRPbdjVRWdtEbWNrh98fH2f07pVI79QEslITKcxOJSsli8zUBDJTEshISSAjOZ7MlATSkwP3M1ISSEsKFLYKOHJ0Wu7OuYs7ed4B3w9bIhHZpx0NLWyoqqe8up7y6gbKqxuoqGmgorqBTTsaqK5v+cr3pCb6yM1IJictiZGHZXDK0CT69kokOz2J7LQk+vRKpG+vRPqkJZKeFK+CjhGeXRVSRDq2o6GFNdvqWLNtF2sq61i7vZ512wNfdzR8ubzTkuLJy0phQFYyYwuyGJCZTP/MFA7LTCY3I5ncjCTSkxM8+i8RL6ncRTxSU9/Myi27KNlSy8rNtZRu3UVp5S4qa5v2jIkzyOudQmHfXpw7uj8FfVIp6JNKfp9UBvZOITMlQTNt6ZDKXaSLOefYUNXAko07WFKxg+WbdrJicy2bdjTuGZOeFM8R/dKYOCyHIf3SGJzdi8E5aRT0SSUxXpeAkgOnchcJs607G1m4oYbPN9SwuHwHi8tr2Bl8AzM+zhjSL43jB/VhZP8Mhh+WzrDcdPpnJmsGLmGlchc5BH6/o2RLLfPXVlG8tpoF66qpqGkAAkU+on86544ZwNF5mRw1IJNhh6WRFO/zOLX0BCp3kQPg9zuWb97Jx6u380nZduatqdozK8/NSOLYw3tzxYRCxhX05sgBGSQnqMjFGyp3kU5s3dnIuysreX/VNj4o3UZVXTMAhX1TOfvo/hxX2Ifxg/owsHeKDq1IxFC5i+zF73csKq9h9vKtzCnZytKNOwHITkti4rAcJgzJ5qQj+jIgK8XjpCL7pnIXAZpb/Xy0ehuzlm7h7eVbqKxtwhdnHFvQm1snj2DisBxG9k/XzFyihspdeqyWNj8flm7jtcWbeGvZFnY0tNAr0cek4f04Y1Qupw7vR2aqPgAk0UnlLj2Kc46FG2p4dWEF/1y8ie11zaQnxXPGqFzOOro/XxuarTdBJSao3KVH2LKzkb99Vs5LC8opq6wjKT6O00fmcv7YAUwanqPTEyXmqNwlZrX5He+tquT5T9cze8VW2vyO4wp7c+0pgzn76P665orENJW7xJwd9S28ULyeZz9eR3l1A317JXL11wbz78flMyi7l9fxRLqFyl1ixurKXfzxgzW88lkFDS1tjB/Uh9vOGsGZow7T9Vmkx1G5S9RbsK6ax95dzVvLt5Dgi+ObY/O47KRCRg3I8DqaiGdU7hKVnHN8vHo7D81exSdlVWSmJHDjqUP47omF5KQneR1PxHMqd4kqzjk+LN3Ob99eyYJ11fRLT+K/zx3FxePzSU3UX2eR3fTbIFFjwboq7p9VwidlVQzITOaeKUdyUVG+zksX6YDKXSJe6dZd/Or1Fby9fAvZaYnccd4oLjm+QOemi+yHyl0i1vZdTTz49iqen7eelAQfPzlzGFeePEiHX0RCoN8SiTitbX7+9Mk6/vetldQ1t3HJ+AJuOn0o2Wl6o1QkVCp3iSjz1lRx+9+XsGJzLScPyeaO80YxNDfd61giUUflLhGhpr6Ze2eu4IXiDeRlpfDopcfyjSNzdYldkYOkchdPOef4x+JN3DVjKTUNLVw7cTA3nTZUx9VFDpF+g8QzlbVN/OLVL5i1dAtjBmby3FXH61OlImGichdPvLZ4I//96hLqmtu47awRfO/kQcT7dP0XkXBRuUu3qm1s4Y6/L+XlhRWMyc/iNxeNZkg/vWEqEm4qd+k2C9ZVcfMLi6iobuCm04byg68P0WxdpIuo3KXL+f2Ox98v4/5ZJQzISuav153IsYf38TqWSExTuUuXqq5r5pa/fs7sFVs55+j+3PtvR5OhFZBEulxI/yY2s8lmVmJmpWZ2WwfPF5jZHDNbaGaLzezs8EeVaLOkYgfn/v4DPli1jXumHMnDl4xTsYt0k05n7mbmA6YBZwDlwHwzm+GcW9Zu2C+AF51zj5jZKGAmUNgFeSVKvPxZOT97+Quy05J46foTGT0wy+tIIj1KKIdlxgOlzrkyADObDkwB2pe7A3afoJwJbAxnSIkerW1+fjlzOU99uJYTB/fl4UvG0VfXhBHpdqGUex6wod3jcuD4vcbcCbxpZj8AegGnd/SDzOwa4BqAgoKCA80qEW5nYws/eH4h766s5MoJg/j52SN0NoyIR0L5zevo4h5ur8cXA0875wYCZwPPmdlXfrZz7nHnXJFzrignJ+fA00rE2lBVz4WPfMSHpdu491tHc/t5o1TsIh4KZeZeDuS3ezyQrx52uQqYDOCc+9jMkoFsYGs4QkpkW1xew5VPz6ep1c8zV45nwpBsryOJ9HihTK3mA0PNbJCZJQJTgRl7jVkPnAZgZiOBZKAynEElMr27spKpj39CUryPV244ScUuEiE6nbk751rN7EZgFuADnnTOLTWzu4Fi59wM4BbgCTP7EYFDNpc75/Y+dCMx5pWF5fz0r4sZmpvO01ccR25GsteRRCQopA8xOedmEji9sf2229vdXwZMCG80iWRPf7iGO/+xjBMH9+Wx7x6r89dFIow+oSoHbNqcUu6fVcKZo3L5/SXjtFC1SARSuUvInHPcP6uEP8xdzQVjB3D/RWNI0BkxIhFJ5S4hcc7xq9dX8Nh7ZVw8voBfXnAUcXFaAk8kUqncpVPOOX71RqDY/+OEw7l7ypFa21Qkwunf1LJfzjl+/UYJj71bxqUnFKjYRaKEyl3263fvrOLRd1fzneMLuPv8o1TsIlFC5S779McP1vDg26u48NiB3DNFx9hFoonKXTr04vwN3PPaMs466jB+9a2jVewiUUblLl/xxpLN3PbyYk4ZlsODU8fqAmAiUUi/tfIl89dW8cPpCxmTn8Vjlx6rDyiJRCmVu+yxakst33ummIG9U3jysuNISVSxi0QrlbsAsGVnI5c9OY/E+DieuWI8vXsleh1JRA6Byl2ob27lqmfmU9PQwlOXH0d+n1SvI4nIIVK593B+v+PHL3zOso07+f3F4zgqL9PrSCISBir3Hu6BN0t4Y+lm/uucUZw2MtfrOCISJir3Huzlz8r5w9zVXHJ8AVdOKPQ6joiEkcq9h1pcXsNtL3/BiYP7ctf5ul6MSKxRufdA23Y1ce1zC8hJS+LhS8bpmuwiMUiX/O1hWtr83PDnz6iub+al606ib1qS15FEpAuo3HuY/5m5nHlrqvjd1LE6M0Ykhunf4z3IzC828dSHa7n8pEKmjM3zOo6IdCGVew+xZlsd//nSYsbmZ/Hzs0d6HUdEupjKvQdobGnj+j8tIN5nTPvOMSTG649dJNbpmHsPcNc/lrFicy1PXXEceVkpXscRkW6gKVyMm/nFJv4ybz3XTTyCU4f38zqOiHQTlXsMK6+u57a/LWZMfha3nDnM6zgi0o1U7jGqtc3PzdMX4Xfw0NSx+qCSSA+jY+4x6uE5pRSvq+bBfx/L4X17eR1HRLqZpnMxaNGGGn4/u5QLxg7ggnE6n12kJ1K5x5iG5jZ+/MIi+qUncdeUo7yOIyIeCanczWyymZWYWamZ3baPMd82s2VmttTMng9vTAnVr99YQdm2Oh64aAyZKQlexxERj3R6zN3MfMA04AygHJhvZjOcc8vajRkK/AyY4JyrNjOdc+eB91dV8vRHa7liQiEThmR7HUdEPBTKzH08UOqcK3PONQPTgSl7jbkamOacqwZwzm0Nb0zpTG1jC7e+tJgjcnpx6+QRXscREY+FUu55wIZ2j8uD29obBgwzsw/N7BMzmxyugBKa/5m5gs07G/nNt8eSnODzOo6IeCyUUyE7WqLHdfBzhgKTgIHA+2Z2lHOu5ks/yOwa4BqAgoKCAw4rHfuwdBt/mbeea08ZzNj8LK/jiEgECGXmXg7kt3s8ENjYwZi/O+danHNrgBICZf8lzrnHnXNFzrminJycg80s7dQ1tXLr3xYzOLsXPzpDn0IVkYBQyn0+MNTMBplZIjAVmLHXmFeBUwHMLJvAYZqycAaVjt0/q4SKmgbuu3C0DseIyB6dlrtzrhW4EZgFLAdedM4tNbO7zez84LBZwHYzWwbMAX7qnNveVaElYMG6ap75eC2XnVhIUWEfr+OISAQx5/Y+fN49ioqKXHFxsSevHQta2vyc+9AH7Gxs4a0fTyQtSVeSEOkJzGyBc66os3FqhCj1xPtllGyp5YnvFqnYReQrdPmBKLRuex2/e3sVZx11GGeMyvU6johEIJV7lHHO8YtXl5Doi+PO84/0Oo6IRCiVe5T55xebeH/VNn46eTi5GclexxGRCKVyjyK7mlq557VlHJWXwXeOP9zrOCISwfROXBR56J1VbNnZxKOXHosvrqMPDouIBGjmHiVWbqnlyQ/WMPW4fMYV9PY6johEOJV7FHDO8d+vLiEtOZ7/1BUfRSQEKvco8M8vNvHpmip++o3h9OmV6HUcEYkCKvcI19Dcxr0zVzCqfwZTj9OVNEUkNCr3CPfYe6upqGngjvNG6U1UEQmZyj2CVdQ08Oi7qzlndH+OH9zX6zgiEkVU7hHs3pnLcQ5+fvZIr6OISJRRuUeo4rVVvLZ4E9dOPIK8rBSv44hIlFG5RyDnHP/vn8vJzUjiuomDvY4jIlFI5R6B/rF4E4s21HDLmcNJTdSHiEXkwKncI0xjSxu/fn0FI/tn8G/HDPQ6johEKZV7hHn6o7VU1DTwi3NG6tRHETloKvcIUlXXzLTZpZw6PIcJQ7K9jiMiUUzlHkGmzSmlrrmVn+nURxE5RCr3CLGhqp7nPl7HhccOZFhuutdxRCTKqdwjxG/fWokZ3Hz6MK+jiEgMULlHgOWbdvLKogoun1DIAH1gSUTCQOUeAe57YwXpSfHcMHGI11FEJEao3D32adl25pRUcsOpQ8hMTfA6jojECJW7h5xz3D+rhNyMJC4/qdDrOCISQ1TuHpq7spLiddX84OtDSU7weR1HRGKIyt0jfr/jgVklFPRJ5dtF+V7HEZEYo3L3yBtLN7N0405uPn0oifH6YxCR8FKreKDN7/jNmyUM7ZfGlLF5XscRkRikcvfAqwsrWF1Zxy1nDtPFwUSkS4RU7mY22cxKzKzUzG7bz7gLzcyZWVH4IsaW1jY/D81exZEDMvjGkYd5HUdEYlSn5W5mPmAacBYwCrjYzEZ1MC4d+CHwabhDxpKXF1awbns9Pzp9GGaatYtI1whl5j4eKHXOlTnnmoHpwJQOxt0D3Ac0hjFfTGlp8/PQO6sYPTCT00b28zqOiMSwUMo9D9jQ7nF5cNseZjYOyHfOvRbGbDHnpQXllFc3aNYuIl0ulHLvqIXcnifN4oDfArd0+oPMrjGzYjMrrqysDD1lDGhu9fPw7FLG5mcxaXiO13FEJMaFUu7lQPtP2QwENrZ7nA4cBcw1s7XACcCMjt5Udc497pwrcs4V5eT0rIL764INVNQ08KMzNGsXka4XSrnPB4aa2SAzSwSmAjN2P+mc2+Gcy3bOFTrnCoFPgPOdc8VdkjgKNbf6+cOc1YwryOKUoVo+T0S6Xqfl7pxrBW4EZgHLgRedc0vN7G4zO7+rA8aClz8rp6KmgZtOG6pZu4h0i/hQBjnnZgIz99p2+z7GTjr0WLGjpc3Pw3NKGTMwk4nDetahKBHxjj6h2sVe+ayC8uoGbjpds3YR6T4q9y7UGpy1H52XyanDdV67iHQflXsX+vuijayvqueHOtYuIt1M5d5F2vyOaXNLGdk/g9P1aVQR6WYq9y7y+pJNlFXW8YOvD9GsXUS6ncq9C/j9jodnlzKkXxqTdeVHEfGAyr0LvLNiKys21/L9U48gTtdrFxEPqNzDzDnHw7NXUdAnlfNGD/A6joj0UCr3MHt/1TY+L9/B9ZOOIN6n3Ssi3lD7hNnDc0o5LCOZbx2jtVFFxDsq9zCav7aKeWuquOaUwSTF+7yOIyI9mMo9jKbNKaVPr0QuHl/gdRQR6eFU7mGypGIHc0squerkQaQkatYuIt5SuYfJI3NXk54Uz6UnHO51FBERlXs4rK7cxcwlm/iPEw8nMyXB6zgiIir3cHh07mqS4uO48uRBXkcREQFU7oesoqaBVxZWMPW4ArLTkryOIyICqNwP2RPvlQFw9SmDPU4iIvIvKvdDsH1XE9Pnr2fK2DzyslK8jiMisofK/RA8/dFamlr9XD9Js3YRiSwq94NU29jCMx+t5cxRuQzpl+51HBGRL1G5H6TnP13PzsZWbpg0xOsoIiJfoXI/CI0tbfzfB2s4eUg2Y/KzvI4jIvIVKveD8LfPyqmsbeKGSUd4HUVEpEMq9wPU2ubnsXfLGJOfxYlH9PU6johIh1TuB+ifX2xifVU9N0w6Qgtfi0jEUrkfAOccj8xdzdB+aZwxMtfrOCIi+6RyPwBzSgILX183UQtfi0hkU7mHyDnHH+asJi8rhfPHauFrEYlsKvcQzVtTRfG6aq6dOJgELXwtIhFOLRWiaXNXk52WyLeL8r2OIiLSqZDK3cwmm1mJmZWa2W0dPP9jM1tmZovN7B0zi6nliL4o38F7Kyu56uTBJCdoCT0RiXydlruZ+YBpwFnAKOBiMxu117CFQJFzbjTwEnBfuIN66Q9zS0lPjufSE7TwtYhEh1Bm7uOBUudcmXOuGZgOTGk/wDk3xzlXH3z4CTAwvDG9U7q1ljeWbubykwpJT9YSeiISHUIp9zxgQ7vH5cFt+3IV8HpHT5jZNWZWbGbFlZWVoaf00CNzy0iO93HFBC2hJyLRI5Ry7+iEbtfhQLNLgSLg/o6ed8497pwrcs4V5eTkhJ7SIxuq6nl1UQVTx+fTp1ei13FEREIWH8KYcqD9KSIDgY17DzKz04H/AiY655rCE89bj7y7Gp8Z156iC4SJSHQJZeY+HxhqZoPMLBGYCsxoP8DMxgGPAec757aGP2b327SjgZeKy7moaCCHZSZ7HUdE5IB0Wu7OuVbgRmAWsBx40Tm31MzuNrPzg8PuB9KAv5rZIjObsY8fFzUef6+MNue4bqJm7SISfUI5LINzbiYwc69tt7e7f3qYc3mqsraJv8xbzzfH5ZHfJ9XrOCIiB0yfUO3A/31QRnOrX4txiEjUUrnvpaqumT99vI5zRg9gcE6a13FERA6Kyn0vT7xfRn1LGz/8uha+FpHopXJvp6qumWc+Wst5owcwNDfd6zgiIgdN5d7OE++X0dDSxg9P06xdRKKbyj1o+66mPbP2If00axeR6KZyD3ri/TXBWftQr6OIiBwylTuB89qf/Xgt548ZwJB+OkNGRKKfyh2YNqeUplY/N58+zOsoIiJh0ePLfUNVPX/+dB3fLspnUHYvr+OIiIRFjy/3B99ehZnpDBkRiSk9utxXbanllYXlXHbi4fTPTPE6johI2PTocn/gzRJSE+O5fpJm7SISW3psuc9fW8WspVu4+muDtcqSiMScHlnufr/jnteWkZuRxNWnaG1UEYk9PbLcX11UweLyHdw6eQSpiSFd0l5EJKr0uHKvb27lvjdKGD0wkwvG5nkdR0SkS/S4cn/s3TI272zk9nNHERdnXscREekSParcy6vreey91Zwzuj9FhX28jiMi0mV6TLk75/jFq0uIM+PnZ4/0Oo6ISJfqMeX+2uJNzC2p5CdnDicvSx9YEpHY1iPKfUd9C3f9YymjB2Zy2UmFXscREelyPeI8wHtfX051fQvPXDken95EFZEeIOZn7u+urGT6/A187+RBHDkg0+s4IiLdIqbLfcvORn78wiJGHJbOj87QtdpFpOeI2XJv8ztumr6Q+uY2Hr5kHMkJPq8jiYh0m5g95v772av4pKyKBy4aowWvRaTHicmZ+9vLtvDQO6v41rg8Ljx2oNdxRES6XcyV+4el27jh+c84emAW91xwlNdxREQ8EVPlvmBdNVc/W8ygvr145orj6JUUs0edRET2K6RyN7PJZlZiZqVmdlsHzyeZ2QvB5z81s8JwB+3MB6u2ccVT8+iXnsRz3xtPVqoW4BCRnqvTcjczHzANOAsYBVxsZqP2GnYVUO2cGwL8Fvh1uIPuS2NLG3fOWMqlf/yUnPQk/vS94+mXntxdLy8iEpFCOW4xHih1zpUBmNl0YAqwrN2YKcCdwfsvAQ+bmTnnXBiz7uGcY9XWXXxStp1nP15H6dZdXH5SIbdOHkFKok55FBEJpdzzgA3tHpcDx+9rjHOu1cx2AH2BbeEI2d70eet54M0Stu1qBqCwbyrPXjmeU4blhPulRESiVijl3tHFWPaekYcyBjO7BrgGoKCgIISX/qrcjGS+NjSHEwb34cTB2eT3ScFM14sREWkvlHIvB/LbPR4IbNzHmHIziwcygaq9f5Bz7nHgcYCioqKDOmRz6oh+nDqi38F8q4hIjxHK2TLzgaFmNsjMEoGpwIy9xswALgvevxCY3VXH20VEpHOdztyDx9BvBGYBPuBJ59xSM7sbKHbOzQD+CDxnZqUEZuxTuzK0iIjsX0if8nHOzQRm7rXt9nb3G4GLwhtNREQOVkx9QlVERAJU7iIiMUjlLiISg1TuIiIxSOUuIhKDzKvT0c2sElh3kN+eTRdc2qALRVPeaMoK0ZU3mrJCdOWNpqxwaHkPd851er0Vz8r9UJhZsXOuyOscoYqmvNGUFaIrbzRlhejKG01ZoXvy6rCMiEgMUrmLiMSgaC33x70OcICiKW80ZYXoyhtNWSG68kZTVuiGvFF5zF1ERPYvWmfuIiKyH1FX7p0t1h1JzGytmX1hZovMrNjrPHszsyfNbKuZLWm3rY+ZvWVmq4Jfe3uZcbd9ZL3TzCqC+3eRmZ3tZcb2zCzfzOaY2XIzW2pmNwW3R9z+3U/WiNy/ZpZsZvPM7PNg3ruC2weZ2afBfftC8BLlkZr1aTNb027fjg37izvnouZG4JLDq4HBQCLwOTDK61z7ybsWyPY6x37ynQIcAyxpt+0+4Lbg/duAX3udcz9Z7wR+4nW2feTtDxwTvJ8OrCSwwHzE7d/9ZI3I/Utg5be04P0E4FPgBOBFYGpw+6PA9RGc9Wngwq587Wibue9ZrNs51wzsXqxbDoJz7j2+umLWFOCZ4P1ngAu6NdQ+7CNrxHLObXLOfRa8XwssJ7DWcMTt3/1kjUguYFfwYULw5oCvAy8Ft0fKvt1X1i4XbeXe0WLdEfuXkMAf4ptmtiC4fmw0yHXObYLALz0Q6Wsa3mhmi4OHbTw/xNERMysExhGYtUX0/t0rK0To/jUzn5ktArYCbxH4F32Nc641OCRiumHvrM653fv2l8F9+1szSwr360ZbuYe0EHcEmeCcOwY4C/i+mZ3idaAY8whwBDAW2AT8xts4X2VmacDfgJudczu9zrM/HWSN2P3rnGtzzo0lsKbzeGBkR8O6N1XH9s5qZkcBPwNGAMcBfYBbw/260VbuoSzWHTGccxuDX7cCrxD4SxjptphZf4Dg160e59kn59yW4C+OH3iCCNu/ZpZAoCz/7Jx7Obg5IvdvR1kjff8COOdqgLkEjmNnmdnu1eUirhvaZZ0cPBTmnHNNwFN0wb6NtnIPZbHuiGBmvcwsffd94Exgyf6/KyK0X+z8MuDvHmbZr90lGfRNImj/mpkRWFt4uXPuf9s9FXH7d19ZI3X/mlmOmWUF76cApxN4n2AOcGFwWKTs246yrmj3P3gj8N5A2Pdt1H2IKXg61oP8a7HuX3ocqUNmNpjAbB0Ca9U+H2lZzewvwCQCV6jbAtwBvErgrIMCYD1wkXPO8zcy95F1EoFDBo7AmUnX7j6e7TUzOxl4H/gC8Ac3/5zAseyI2r/7yXoxEbh/zWw0gTdMfQQmqC865+4O/s5NJ3BDV9uOAAAAVklEQVSYYyFwaXBm7Jn9ZJ0N5BA41LwIuK7dG6/hee1oK3cREelctB2WERGREKjcRURikMpdRCQGqdxFRGKQyl1EJAap3EVEYpDKXUQkBqncRURi0P8HXlB5ntqUsjMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(*G.step())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to find out what manipulations must be made (changes to $u$) in order to get the system to follow a specific desired trajectory (which we will call $r$ for the reference trajectory). We will allow the controller to make a certain number of moves. This is called the control horizon, $M$. We will the observe the effect of this set of moves (called a \"move plan\") for time called the prediction horizon ($P$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Controller parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 10  # Control horizon\n",
    "P = 20  # Prediction horizon\n",
    "DeltaT = 1  # Sampling rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "tcontinuous = numpy.linspace(0, P*DeltaT, 1000)  # some closely spaced time points\n",
    "tpredict = numpy.arange(0, P*DeltaT, DeltaT)   # discrete points at prediction horizon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose a first order setpoint response similar to DS or Dahlin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau_c = 1\n",
    "r = 1 - numpy.exp(-tpredict/tau_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an initial guess we choose a step in $u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = numpy.ones(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initital state is zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = numpy.zeros(G.to_ss().A.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extend(u):\n",
    "    \"\"\"We optimise the first M values of u but we need P values for prediction\"\"\"\n",
    "    return numpy.concatenate([u, numpy.repeat(u[-1], P-M)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prediction(u, t=tpredict, x0=x0):\n",
    "    \"\"\"Predict the effect of an input signal\"\"\"\n",
    "    t, y, x = scipy.signal.lsim(G, u, t, X0=x0, interp=False)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1475ccf8>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl81NW9//HXh30LCZAEEBIIO8hOCCC2RUFF6lKxKmgVAcW6VWtrL79rXXu9127a25arUmUREVcsWLFKryJWBBPCIoQthABhS8KeQMgy5/dHBm4aEzJAJt9k5v18POaRmfmezHweX77z5psz53uOOecQEZHQUs/rAkREpPop3EVEQpDCXUQkBCncRURCkMJdRCQEKdxFREKQwl1EJAQp3EVEQpDCXUQkBDXw6o2jo6Nd586dvXp7EZE6afXq1bnOuZiq2nkW7p07dyYlJcWrtxcRqZPMbGcg7dQtIyISghTuIiIhSOEuIhKCFO4iIiFI4S4iEoIU7iIiIUjhLiISgjwb5y4iEi6OFxSRmXuCjNw8duTmM7pXW/p1jAzqeyrcRUSqQUFRCbsOnWBHbn7pLaf0Z0ZuPrl5p860M4PoFo0V7iIitYXP59hz5CTbc0rPwDP94b0jN589R07i3P+1jW7RmIToZlzeK4aE6BYkRDenS0xz4ls3o0nD+kGvVeEuIlKBoyeK2Lz/GJv3H/ffjrF1/3HyC0vOtIlo3ICEmOYM6dSKGwd3pEtMcxKim9M5ujktmzT0sHqFu4iEucJiHxm5eWzZf5xN+0pDfMv+4+w7WnCmTVSzhvRsG8FNiXH0bBdBt9gWdG7TnOgWjTAzD6uvnMJdRMLG4fxC1mUdYdO+42zxn5Vvz8mjqKS0P6VhfaNrTAuGd2lDz3YR9GoXQa92LWnbsnGtDfHKKNxFJCSdKi5h495jrNt9hLX+286DJ85svyiyCT3bRXBZr9gzId4lpjkN64fGCHGFu4jUec45duTms3b3kTNhnrbv2Jkz8rYtGzMwLooJQ+MZEBfJxe0jiWzmbZ94sCncRaTOOZRfyNrdh1m7++iZQD96sgiAZo3q069DJFMuTWBQXBQD41rRLrKJxxXXPIW7iNR62ccK+CrjICvSD7Jyx8Ez3Sv1DHq0jeDqvu0YGBfFwPgousdGUL9e3eofDwaFu4jUOofzC1mZcZAV2w+yYnsu23PyAWjZpAHDurRhYlI8A+Oi6NchkuaNFWMV0V4REc8dKygiecchf5gfZNO+Y0BpF0tSQmtuGRrHiC7R9Lmopc7KA6RwF5Ead7KwhJSd/xfm32QdweegUYN6JHZqxc+v7MGIrm3o3zEqZEav1DSFu4jUiO05efwj7QD/uzmbNbsOU1TiaFDPGBgXxf2XdWNE1zYMjm9VI5fmhwOFu4gERYnPsWbXYZZuOsDStANk+PvNe7dvyZSRCYzo2oahnVurzzxItFdFpNqcLCzhi205LE07wKebszmYX0iDesbwLm2YNKIzo3vH0rFVM6/LDAsKdxG5IDnHT/Hp5tKz8y+25XKq2EdEkwZc1jOWMX3aMqpnjOeTaIUjhbuInBPnHNtz8lials3StP2s2X0E56BDVFMmJsUzpndbkhJa06iBvgj1ksJdRAKy6+AJFq7JYvHavWTklvaf9+3QkodH9+CKPm3p3T6izk2uFcoU7iJSqaMni/hw/T7eX5NFcuZhzGB4Qhsmj+zMmD5taR/Z1OsSpRIKdxH5F0UlPj7fksP7a/awdNMBCot9dIttwS/G9uQHAztwUZQCvS5QuIsIzjm+2XOUhal7+GDdXg7mF9K6eSNuTYpn/OAO9OsQqS6XOkbhLhLG9h45yV/X7mFh6h7Ss/NoVL8eY/rEMn5QR77XM0ZXh9ZhCneRMJN3qpi/b9jPwtQsvso4iHOQ2KkV/3lDP77fr33Iz3MeLhTuImEiPTuPOSt2sDB1DycKS4hv3YyHRnfnhkEd6NSmudflSTVTuIuEMJ/P8fm2HGZ/mcnyrTk0alCP6wZcxIShcQzp1Er96CFM4S4SgvJPFbMwNYvZKzLJyMknJqIxj1zRg1uHxRPdorHX5UkNULiLhJDdh07w2leZvJm8m+MFxQzoGMkfbhnIuH7tdcVomAko3M1sLPDfQH3gFefcc+W2xwNzgSh/m+nOuSXVXKuIVMA5x9c7DjH7y0w+SduPmTG2bzumjExgcHyUul7CVJXhbmb1gRnAFUAWkGxmi51zaWWa/RJ42zn3opn1AZYAnYNQr4j4FRSV8MG6vcz+MpO0fceIataQe77XlduHd9KFRhLQmXsSkO6cywAwszeB64Gy4e6Alv77kcDe6ixSRP5P9vECXl+5izdW7SQ3r5AebVvwX+P78YOBHWjaSAtdSKlAwr0DsLvM4yxgWLk2TwGfmNmDQHNgTLVUJyJn7D9awIzP0nkreTdFPh+X94xlyqUJXNK1jbpe5FsCCfeKjhpX7vFEYI5z7vdmNgKYZ2Z9nXO+f3khs2nANID4+PjzqVck7OTmneLFZduZt3InPp/j5qFxTPtOFzpHa2y6VC6QcM8C4so87si3u12mAmMBnHNfmVkTIBrILtvIOTcTmAmQmJhY/j8IESnjcH4hLy/PYO6KTE4Vl3Dj4I78ZHR34lprJSOpWiDhngx0N7MEYA8wAbi1XJtdwGhgjpn1BpoAOdVZqEi4OHqyiFf/uYNZ/9xBfmEx1w24iIdGd6dLTAuvS5M6pMpwd84Vm9kDwMeUDnOc5ZzbaGbPACnOucXAz4C/mNlPKe2yudM5pzNzkXOQd6qYOV/uYObyDI4VFDOuXzseHtODHm0jvC5N6qCAxrn7x6wvKffcE2XupwEjq7c0kfBwsrCE177K5OXlGRzKL2RM71h+ekUPLr4o0uvSpA7TFaoiHikoKmHB17uY8dl2cvNO8d0eMTxyRQ8GxkV5XZqEAIW7SA0rLPbxzurd/PnTdPYdLWB4l9a8+KPBDO3c2uvSJIQo3EVqiHOOjzfu59klm9h96CRDOrXi9zcN4JJu0V6XJiFI4S5SA3bk5vPk4o0s35pDr3YRzJ48lFE9YnTxkQSNwl0kiE4WljDjs3RmLs+gcYN6PHltH24f3okGWr5OgkzhLhIEzjk+STvAMx+ksefIScYP6sD0cb2IjWjidWkSJhTuItUsMzefpz7YyLItOfRsG8Fb04YzrEsbr8uSMKNwF6kmJwtLeHFZOi99nkGjBvV4/Jo+3DGiEw3VBSMeULiLVIOlaQd4+oONZB0+yfUDL+Kxcb2JbakuGPGOwl3kAuw8mM/TH6Tx6eZsuse2YMHdwxnRVV0w4j2Fu8h5KCgq4cVl23nx8+00rGc8Nq43d47srC4YqTUU7iLn6NPNB3hqcRq7Dp3g2gGlXTDtItUFI7WLwl0kQMcKinhy0UbeX7OHbrEteOOuYbq6VGothbtIAFZmHORnb69j/7ECHhrdnfsv60ajBuqCkdpL4S5yFqeKS3h+6VZmLs+gU+tmvPvjEQyKb+V1WSJVUriLVGLrgeM89OZaNu07xsSkeH75/d40b6yPjNQNOlJFyvH5HHNWZPLc3zcT0bgBr9yRyJg+bb0uS+ScKNxFyth/tIBH313HF9tyGd0rludu7E9MRGOvyxI5Zwp3Eb8P1+/j39//hsJiH8/e0Jdbk+I1Ja/UWQp3CXvHCop4atFGFq7Zw4C4KF64eQBdYlp4XZbIBVG4S1j7eschfvrWWvYfK+Ano7vz4OXddJWphASFu4SlwmIfzy/dysvLtxPfuhnv/HgEgzXEUUKIwl3Czjb/EMe0fceYMDSOx6/poyGOEnJ0REtYeXd1Fo+9/w3NGzfgL3ckcoWGOEqIUrhLWCgq8fHsh5uYsyKTS7q24Q8TBmrJOwlpCncJeYfyC7l/fipfZRxkysgE/n1cLy1QLSFP4S4hbePeo0x7bTU5eaf4/U0DuHFIR69LEqkRCncJWR+s28uj764jqmkj3rlnBAPiorwuSaTGKNwl5JT4HL/9eAsvfb6dxE6tePFHQzSFgIQdhbuElKMni3jozTUs25LDrcPieeraizXvuoQlhbuEjG0HjjNt3mqyDp/g2Rv6ctuwTl6XJOIZhbuEhE827ueRt9fRpGF93rh7OEM7t/a6JBFPKdylTvP5HH/6NJ0X/rGV/h0jefn2IbSPbOp1WSKeC6gz0szGmtkWM0s3s+mVtLnZzNLMbKOZvVG9ZYp8W96pYu6dv5oX/rGV8YM68PY9IxTsIn5VnrmbWX1gBnAFkAUkm9li51xamTbdgf8HjHTOHTaz2GAVLAKQmZvPtHkpbM/J5/Fr+jBlZGfNvS5SRiDdMklAunMuA8DM3gSuB9LKtLkbmOGcOwzgnMuu7kJFTvt8aw4PvpFKvXrGa1OSGNkt2uuSRGqdQLplOgC7yzzO8j9XVg+gh5l9aWYrzWxsdRUoUtbrK3cyefbXXBTVlMX3X6pgF6lEIGfuFf2t6yp4ne7AKKAj8IWZ9XXOHfmXFzKbBkwDiI+PP+diJXw5V/rF6fNLt3J5r1j+NHGQpukVOYtAztyzgLgyjzsCeytos8g5V+Sc2wFsoTTs/4VzbqZzLtE5lxgTE3O+NUuY8fkcT3+QxvNLtzJ+cAdevn2Igl2kCoGEezLQ3cwSzKwRMAFYXK7NX4HLAMwsmtJumozqLFTCU1GJj0feXsucFZlMvTSB3/1wgJbBEwlAlac/zrliM3sA+BioD8xyzm00s2eAFOfcYv+2K80sDSgBHnXOHQxm4RL6ThaWcO/81SzbksOjV/XkvlFdNSJGJEDmXPnu85qRmJjoUlJSPHlvqf2OnChk6twU1uw6zLM39GNikr6jEQEws9XOucSq2qnjUmqdA8cKuOPVr9mRm8+MWwdzdb/2XpckUuco3KVW2ZGbz+2vruJwfiFzJg/lEg11FDkvCnepNTbsOcqds7/G52DBtOH076jFNUTOl8JdaoWvth/k7tdSiGzakNemJtE1poXXJYnUaQp38dzHG/fz4II1xLduxrypSZr8S6QaKNzFU2+n7Gb6e+vp3zGK2XcOpVXzRl6XJBISFO7imZc/385/fbSZ73SP5qUf6apTkeqkT5PUOOccz320mZeXZ3BN//Y8f/NArXMqUs0U7lKjikt8/Pv73/B2Sha3D+/EU9ddTP16uupUpLop3KXGlPgcP3tnHYvW7uWh0d15eEx3TScgEiQKd6kRPp9j+nvrWbR2L49e1ZP7L+vmdUkiIU0dnRJ0zjkeX7SBd1Zn8dDo7gp2kRqgcJegcq50Lvb5q3Zx76iuPDzmW9P8i0gQKNwlaE6Pijk9F/svruqpPnaRGqJwl6B5YelWXl6ewe3DO/HL7/dWsIvUIIW7BMWfP93GHz9NZ8LQOJ6+7mIFu0gNU7hLtZu5fDu/+2Qr4wd14D9v6Ec9jWMXqXEKd6lWs7/cwX8u2cw1/dvzmx/2V7CLeEThLtVm/qqdPP1BGldd3JYXbhlIAy1kLeIZffqkWryTspvH3t/A5b1i+dPEwTRUsIt4Sp9AuWCL1u7hF++t5zvdo/mf2wZrEjCRWkCfQrkgS77ZxyNvr2N4Qhtm3p5Ik4b1vS5JRFC4ywVYmnaAnyxYw6C4KF6ZlEjTRgp2kdpC4S7nZdmWbO6fn8rFHSKZPXmoFtoQqWUU7nLOvkzPZdq81XRv24LXJicR0aSh1yWJSDkKdzknq3ceYurcZLpEN+f1qcOIbKZgF6mNFO4SsPTsPKbMSaF9ZFNev2uYFrMWqcUU7hKQ7GMFTJr1NQ3rG3MnJxHdorHXJYnIWehbMKnS8YIi7pydzOEThbw1bQTxbZp5XZKIVEHhLmdVWOzjvvmpbDlwnFcnJdKvY6TXJYlIANQtI5VyzjF94Xq+2JbLc+P7MapnrNcliUiAFO5Sqd99soWFqXv42RU9uCkxzutyROQcKNylQvNW7mTGZ9uZmBTPA5drQWuRuiagcDezsWa2xczSzWz6Wdr90MycmSVWX4lS0z7ZuJ8nF21gTO9YfnW9VlESqYuqDHczqw/MAK4G+gATzaxPBe0igJ8Aq6q7SKk5q3ce5sEFa+jfMYo/ThykOdlF6qhAPrlJQLpzLsM5Vwi8CVxfQbtfAb8BCqqxPqlB23PyuGtuMu0jm/DqpESaNdJgKpG6KpBw7wDsLvM4y//cGWY2CIhzzv3tbC9kZtPMLMXMUnJycs65WAme7OOlFynVr2fMnZJEG12kJFKnBRLuFXW4ujMbzeoBLwA/q+qFnHMznXOJzrnEmJiYwKuUoMo7VcyUOckcyi9k1p1D6dSmudclicgFCiTcs4Cy4+A6AnvLPI4A+gLLzCwTGA4s1peqdUNRSelFSpv2HWfGbYPp3zHK65JEpBoEEu7JQHczSzCzRsAEYPHpjc65o865aOdcZ+dcZ2AlcJ1zLiUoFUu1cc4x/b1vWL41h/+6oR+X6SIlkZBRZbg754qBB4CPgU3A2865jWb2jJldF+wCJXieX7qV91KzeHhMd24eqouUREJJQMMhnHNLgCXlnnuikrajLrwsCbb5q3byp0/TmTA0jodGd/e6HBGpZhrEHIaWph3g8b9u4PJesfzHD/rqIiWREKRwDzNrdx/hwQWp9OsQyZ9v1UVKIqFKn+wwsvfISe6am0JMRGNevXOoLlISCWH6dIeJ/FPFTJ2bwqmiEhbcPUwrKYmEOIV7GPD5HD99ay1b9h9j1p1D6d42wuuSRCTI1C0TBn77yRY+STvA49f00YIbImFC4R7i3l2dxYvLtnPbsHjuvKSz1+WISA1RuIew5MxD/L+F6xnZrQ1PXad52UXCicI9RO0+dIJ75q0mrlUz/ufWITTUkEeRsKJPfAg6XlDE1LnJlPgcr0xKJLJZQ69LEpEapnAPMSU+x4ML1pCRk8+Ltw2mS0wLr0sSEQ9oKGSIefbDTSzbksOzN/Tlkm7RXpcjIh7RmXsIeWPVLmZ9uYPJIztz27BOXpcjIh5SuIeIFdtzeWLRBkb1jOGxcb29LkdEPKZwDwEZOXnc+3oqCdHN+eNETQYmIgr3Ou/oiSLumptC/XrGq5OG0rKJRsaIiMK9Tisq8XHfG6vZffgEL/1oCPFtmnldkojUEhotU0c553hq8Ua+TD/Ib3/Yn6SE1l6XJCK1iM7c66i5KzKZv2oX93yvCzclav1TEflXCvc6aNmWbJ75WxpX9GnLv13Vy+tyRKQWUrjXMenZx3nwjTX0bNeSP9wykHr1NBmYiHybwr0OOZxfyNS5KTRuWJ9XJiXSvLG+MhGRiikd6ojCYh/3zl/NvqMFvDltOB2imnpdkojUYjpzrwOcczy5eCMrMw7x6xv7MTi+ldcliUgtp3CvA+asyGTB17u4d1RXbhjU0etyRKQOULjXcp9vzeFX/pExj17Z0+tyRKSOULjXYunZeTzwRio92kZoZIyInBOFey115EQhd81NpnGDehoZIyLnTIlRCxWV+Lhvfip7jxSwYNowOrbSnDEicm4U7rXQ0x9sZMX2g/z+pgEM6aQ5Y0Tk3KlbppZ57atMXl9ZOmfMjUM0MkZEzo/CvRb5YlsOT3+QxuhesfxCc8aIyAUIKNzNbKyZbTGzdDObXsH2R8wszczWm9n/mpkW8DxHGTl53D8/lW4xLfjviYOor5ExInIBqgx3M6sPzACuBvoAE82sT7lma4BE51x/4F3gN9VdaCg7vZpSg/qlI2NaaGSMiFygQM7ck4B051yGc64QeBO4vmwD59xnzrkT/ocrAXUWB6i4xMf9b6Sy+/AJXr59CHGtNTJGRC5cIOHeAdhd5nGW/7nKTAU+upCiwsmv/pbGP9NzefaGfgztrJExIlI9Avn7v6LOX1dhQ7MfAYnA9yrZPg2YBhAfHx9giaHr9ZU7mfvVTu7+TgI3azUlEalGgZy5ZwFlk6cjsLd8IzMbAzwGXOecO1XRCznnZjrnEp1ziTExMedTb8hYkZ7Lk4s3clnPGKZf3dvrckQkxAQS7slAdzNLMLNGwARgcdkGZjYIeJnSYM+u/jJDy47cfO6dn0qX6Ob8USNjRCQIqgx351wx8ADwMbAJeNs5t9HMnjGz6/zNfgu0AN4xs7VmtriSlwt7R08WMXVuMvUMXp00lIgmDb0uSURCUEBj7pxzS4Al5Z57osz9MdVcV0g6VVzCPfNS2HXwBK/fNYz4NhoZIyLBoQHVNcTnczzy1jpWZhziD7cMZHiXNl6XJCIhTNMP1ADnHM/8LY0Pv9nHY+N684NBZxtJKiJy4RTuNeClzzOYsyKTqZcmcPd3u3hdjoiEAYV7kL27Ootf/30z1w64iMfGacijiNQMhXsQLduSzb+9t56R3drwu5v6a5k8EakxCvcgWbf7CPfNT6Vn2whe+tEQGjeo73VJIhJGFO5BsCM3nylzkmnTohFzpmgsu4jUPIV7Ncs+XsAds1bhgLmTk4iNaOJ1SSIShhTu1SjvVDGTZyeTe7yQVycl0iWmhdcliUiY0kVM1aSw2MeP561m8/7jvHJHIoPiW3ldkoiEMZ25VwOfz/Hou+v4Z3ouz43vx2W9Yr0uSUTCnMK9Gjz3980sWruXR6/qyU2al11EagGF+wV65YsMZi7P4I4RnbhvVFevyxERARTuF2TR2j38x4ebuLpvO5689mLMdJGSiNQOCvfz9M9tufz8nXUkJbTmhVsGasENEalVFO7nYcOeo9wzL4Uu0S34yx2JNGmoq09FpHZRuJ+jHbn53Dk7mcimDZk7JYnIprr6VERqH4X7OVi98xA3vriCEp+PuVOSaBepq09FpHZSuAfow/X7mPiXVbRs0oCF942ke9sIr0sSEamUrlCtgnOOl5dn8NxHm0ns1IqZdyTSunkjr8sSETkrhftZFJf4eHzRRhZ8vYtrB1zEb3/YX1+eikidoHCvxPGCIu5/Yw3Lt+Zw36iu/PzKnlpsQ0TqDIV7BfYdPcnk2clsy87jufH9mJAU73VJIiLnROFezsa9R5kyJ5n8UyXMunMo3+sR43VJIiLnTOFexmebs3ngjVRaNm3IOz8eQe/2Lb0uSUTkvCjc/eat3MmTizbQu31LZt05lLYtNYZdROqusA93n8/x3N83M3N5Bpf3iuVPEwfRvHHY7xYRqePCOsUKikr46Vtr+WjDfm4f3oknr+1Dg/q6rktE6r6wDffcvFPc/VoKa3cf4Zff783USxM0Za+IhIywDPf07Dwmz/ma7GOnePG2wYzt297rkkREqlVYhXtu3ik+WLeXP/xjGw3qGW9OG66FrEUkJIV8uBcUlfCPTQdYmLqHz7fmUOJzDIyL4o8TBhHfppnX5YmIBEVIhrvP50jOPMT7a/bw4fp9HD9VTLuWTbjrOwmMH9SRnu00o6OIhLaAwt3MxgL/DdQHXnHOPVdue2PgNWAIcBC4xTmXWb2lVi0jJ4/31+xhYeoe9hw5SbNG9Rnbtx3jB3VkRNc2WgpPRMJGleFuZvWBGcAVQBaQbGaLnXNpZZpNBQ4757qZ2QTg18AtwSi4vEP5hfxt/V7eS93Dut1HqGcwsls0j17VkysvbkuzRiH5x4mIyFkFknxJQLpzLgPAzN4ErgfKhvv1wFP+++8CfzYzc865aqz1jFPFJXy6KZv3UvewbEs2xT5Hr3YRPDauN9cNvEhXl4pI2Ask3DsAu8s8zgKGVdbGOVdsZkeBNkBudRRZ1lvJu3j2w00cKygmNqIxUy5N4IZBHTQPjIhIGYGEe0Ud1eXPyANpg5lNA6YBxMef3zS67SKbcnmvWMYP7sjIbtHqRxcRqUAg4Z4FxJV53BHYW0mbLDNrAEQCh8q/kHNuJjATIDEx8by6bL7XI0bT8IqIVCGQiVSSge5mlmBmjYAJwOJybRYDk/z3fwh8Gqz+dhERqVqVZ+7+PvQHgI8pHQo5yzm30cyeAVKcc4uBV4F5ZpZO6Rn7hGAWLSIiZxfQOEHn3BJgSbnnnihzvwC4qXpLExGR86X5bUVEQpDCXUQkBCncRURCkMJdRCQEKdxFREKQeTUc3cxygJ3n+evRBGFqg2qk+i6M6rtwtb1G1Xf+OjnnqryS07NwvxBmluKcS/S6jsqovguj+i5cba9R9QWfumVEREKQwl1EJATV1XCf6XUBVVB9F0b1XbjaXqPqC7I62ecuIiJnV1fP3EVE5Cxqdbib2Vgz22Jm6WY2vYLtjc3sLf/2VWbWuQZrizOzz8xsk5ltNLOHKmgzysyOmtla/+2Jil4riDVmmtk3/vdOqWC7mdkf/ftvvZkNrsHaepbZL2vN7JiZPVyuTY3vPzObZWbZZrahzHOtzWypmW3z/2xVye9O8rfZZmaTKmoThNp+a2ab/f9+75tZVCW/e9ZjIcg1PmVme8r8O46r5HfP+nkPYn1vlakt08zWVvK7NbIPq41zrlbeKJ1eeDvQBWgErAP6lGtzH/CS//4E4K0arK89MNh/PwLYWkF9o4C/ebgPM4Hos2wfB3xE6Upaw4FVHv5b76d0/K6n+w/4LjAY2FDmud8A0/33pwO/ruD3WgMZ/p+t/Pdb1UBtVwIN/Pd/XVFtgRwLQa7xKeDnARwDZ/28B6u+ctt/Dzzh5T6srlttPnM/szC3c64QOL0wd1nXA3P9998FRptZjay755zb55xL9d8/DmyidC3ZuuR64DVXaiUQZWbtPahjNLDdOXe+F7VVG+fccr69iljZ42wu8IMKfvUqYKlz7pBz7jCwFBgb7Nqcc58454r9D1dSulKaZyrZf4EI5PN+wc5Wnz87bgYWVPf7eqE2h3tFC3OXD89/WZgbOL0wd43ydwcNAlZVsHmEma0zs4/M7OIaLax0HdtPzGy1f/3a8gLZxzVhApV/oLzcf6e1dc7tg9L/1IHYCtrUhn05hdK/xCpS1bEQbA/4u45mVdKtVRv233eAA865bZVs93ofnpPaHO7VtjB3MJlZC+A94GHn3LFym1Mp7WoYAPwJ+GtN1gaMdM4NBq4G7jez75bbXhv2XyPgOuCdCjZ7vf/Ohaf70sweA4qB+ZU0qepYCKYXga7AQGAfpV0f5Xl+LAITOftZu5f78JzV5nA/l4W5sbMszB0sZtaQ0mCf75xbWH67c+6Ycy7Pf38J0NDMomuqPuelZ7DsAAABzklEQVTcXv/PbOB9Sv/0LSuQfRxsVwOpzrkD5Td4vf/KOHC6u8r/M7uCNp7tS/+Xt9cAtzl/53B5ARwLQeOcO+CcK3HO+YC/VPLenh6L/vwYD7xVWRsv9+H5qM3hXqsX5vb3z70KbHLOPV9Jm3anvwMwsyRK9/fBGqqvuZlFnL5P6RdvG8o1Wwzc4R81Mxw4err7oQZVerbk5f4rp+xxNglYVEGbj4ErzayVv9vhSv9zQWVmY4F/A65zzp2opE0gx0Iwayz7Pc4Nlbx3IJ/3YBoDbHbOZVW00et9eF68/kb3bDdKR3NspfRb9Mf8zz1D6YEM0ITSP+fTga+BLjVY26WU/tm4Hljrv40Dfgz82N/mAWAjpd/8rwQuqcH6uvjfd52/htP7r2x9Bszw799vgMQa/vdtRmlYR5Z5ztP9R+l/NPuAIkrPJqdS+j3O/wLb/D9b+9smAq+U+d0p/mMxHZhcQ7WlU9pXffoYPD167CJgydmOhRrcf/P8x9d6SgO7ffka/Y+/9Xmvifr8z885fdyVaevJPqyum65QFREJQbW5W0ZERM6Twl1EJAQp3EVEQpDCXUQkBCncRURCkMJdRCQEKdxFREKQwl1EJAT9f6Ug0vuWyVYRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tpredict, prediction(extend(u)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(u, x0=x0):\n",
    "    \"\"\"Calculate the sum of the square error for the cotnrol problem\"\"\"\n",
    "    y = prediction(extend(u))\n",
    "    return sum((r - y)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the value of the objective for our step input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.506966650333838"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objective(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we figure out a set of moves which will minimise our objective function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0009187720232727162"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = scipy.optimize.minimize(objective, u)\n",
    "uopt = result.x\n",
    "result.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resample the discrete output to continuous time (effectively work out the 0 order hold value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ucont = extend(uopt)[((tcontinuous-0.01)//DeltaT).astype(int)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the move plan and the output. Notice that we are getting exactly the output we want at the sampling times. At this point we have effectively recovered the Dahlin controller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotoutput(ucont, uopt):\n",
    "    plt.figure()\n",
    "    plt.plot(tcontinuous, ucont)\n",
    "    plt.xlim([0, DeltaT*(P+1)])\n",
    "    plt.figure()\n",
    "    plt.plot(tcontinuous, prediction(ucont, tcontinuous), label='Continuous response')\n",
    "    plt.plot(tpredict, prediction(extend(uopt)), '-o', label='Optimized response')\n",
    "    plt.plot(tpredict, r, label='Set point')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGc9JREFUeJzt3XuspVV5x/Hfc85cuMPgDMPIAAMyWMFWxAkCXmoKwkiMo0YbbKOTakOs2GCiiSCJsW1IsFqa2toqKpE2RNACZaIgIGpJL1wGym2AgQFRRgZmEAWUyzDnPP1jv3vmcNj7XN7nffdaZ+3vJyHnnL3f912LPfv89jprrXctc3cBAMo3kroCAIDBIPABYEgQ+AAwJAh8ABgSBD4ADAkCHwCGBIEPAEOCwAeAIUHgA8CQmJe6AhMtXrzYV6xYkboaADCn3HbbbU+6+5Lpjssq8FesWKH169enrgYAzClm9vOZHEeXDgAMCQIfAIYEgQ8AQ4LAB4AhQeADwJAg8AFgSBD4ADAkspqHv+23L+qC6zbWPv/og/bVqUcf2GCNAKAcWQX+40+/oK/8eJPMZn+uu3TA3gsJfADoI6vAl6RL/vzNessRi2d93rlX3q1rNzzeQo0AoAxF9eG7p64BAOQru8Cv0ZsDAJiB7AK/buLX6fcHgGGSXeBboI1Pjw4A9Jdf4Ndt4dMZBABTCge+mR1sZj8xs/vMbIOZnVU9vr+ZXW9mD1ZfF83oeoG6OKO2ANBXEy38HZI+7e6vk3S8pDPN7ChJZ0u6wd1XSrqh+hkAkEg48N19i7vfXn3/rKT7JB0kaY2ki6vDLpb03plcz2r26ZjRhw8AU2m0D9/MVkh6o6SbJS119y1S50NB0gEzu0aTNQIAdDUW+Ga2l6TLJX3K3Z+ZxXlnmNl6M1sv1e/D53MCAKbWSOCb2Xx1wv4Sd7+ievgJM1tWPb9M0tZe57r7he6+yt1XdY6tXw/GbAGgvyZm6Zikb0m6z90vmPDUOklrq+/XSrpqhleMVgkA0EMTi6e9RdKHJd1tZndUj31O0vmSvmtmH5P0C0kfnMnFas/DN2NaJgBMIRz47v5f6t8sPyl6fQBAM/K70zZ1BQCgUPkFfmDUlg4dAOgvv8BPXQEAKFR+gR9ZHpkmPgD0ld0Wh6lWvXzhpTFddccv9fz2sdrXOOaQRTrm4P0arBUANCe/wE+0PPL/PPSkPnv53aFrHLVsH1191ttC1wCAtmQX+BGRHp3tOzpnX3rG8Xrt0r1nff5nvnenfvHUc4EaAEC7ign8+KJrncDfd/f5WrTnglmfvWBedsMhAPAy2aVUbC2ddKO2LM8MIHf5BX6iQdvuZwXLMwMoVX6BX3vQNm0L28RaPgDyVkzgR3Wjms3QAZQqu8CvK3lXDH34ADKXXeBHWtiRHpVoH75JJD6ArOUX+Klb6gBQqPwCv+55ZvIGmtix8gEgX/kFfrJBW+IaQNmyC/y6beymPidC00KZlgkgY9kFfuxO2zTnStxpCyB/2QU+AKAd2QV+7QZ+sIW969z6XUr06ADIWX6Bz7xMAGhFfoFf+7zYB0V3wLX+FovNTAsFgLbkF/iR3E6Yt3TpAMhdfoFftw+9oZ4gOpQAlCq/wI9My0zcxKeFDyBn2QV+KrsWT6OND6BMxQR+6phmHX0Aucsu8JPdaVt1B9VfPI2lFQDkLcPATztoCwClyi/wA+c20b6eq3vqAsB05qWuwGTJlkdOnNburoe2/U4vjY3XvsaKV+2p3ReMNlgrACXJL/Brr2VjSfvQLTgt85p7HtcnLrk9VIfVRx+or334TaFrAChXdoGfys5pmZEPnECnzq+f2y5JOv/9v6/99pg/6/O/fN0D+s3z22uXD6B82QV+/bVsmq1HKn/0ewfogH12m/V5F/33I8m7pQDkjUHbSedGPnAa2YCFQWMALWkk8M3sIjPbamb3THhsfzO73swerL4umtnFmqjR3LMr7wPTUkl8AFNoqoX/bUmrJz12tqQb3H2lpBuqn6dVvw892sKOpWV4i8No+cExBADlayTw3f1GSU9NeniNpIur7y+W9N6ZXGvu9sVbI33oqbqUAJSvzT78pe6+RZKqrwf0OsjMzjCz9Wa2PlRa8JMi2ocfFezCZxN1ANNKPmjr7he6+yp3XyXN3S78zgdF/ciNrtaZ+j4EAPlrM/CfMLNlklR93TqTk+oHXlDi5ZGTjyEAKF6bgb9O0trq+7WSrprJSdG4TdXKbWqLw9C0VBIfwBSampb5HUn/K+m1ZrbZzD4m6XxJ7zSzByW9s/p5Btdqokaz18jyyKHyd12nXvlGCx/AlBq509bdP9TnqZNme63Ue9qmEl/aQTTxAUwp+aDtKwSDO1XmRQdNo9WmDx/AdPIL/JqiWwzumiXTQGUiIksrkPgAppBd4EcDN1Xmhfvwq7SO9eGT+AD6yy/wE5UbXstGaWfp0MIHMJ38Aj+4p+1cvfkoWm2WVgAwnfwCP1G50T58s2budK1/4xfTMgFMLb/ADwxazmWN3AdAEx/AFPIL/Ohsm9rnxQI3UrbUwF8YgbIBDIfsAn+uim5A0sg8fBr4AKaQXeBH97StG3rRLQabkmoTdQDly24T87nKgoOm8UHjWAv/nl8+rQuuf0Bj4/Uv8idvPkSnHn1g/UoAaFV2gR+ZJRORum0cbZ1Hb/z66cat+vH9W/WG5fvW+ke4f8sz2mvhPAIfyFh+gR8etI3vDVvrvOAsmfA8/OBaPt2G/eV/caLmjc6+p++Uv/9PjTOIAGStmD78sPAm4s38lRCYht9Il9JIzQqMmBH4QObyC/zg+U3csZpSbHnk+uWOB9fyGTHT2Hj98gG0L7vAn6uig6bNLJ4WKH/CdeoYGeHGLyB32QV+dC2dulLPyoz34UfHEDz0GtKlA+Qvv8BPXYGamlqeONkWix577c1MgRmdAAYgv8CvvbRAUxugpFkeOdqlEi/faw/YStKoiRY+kLkMA7+Z4J5r4oPNsb8wxj3WLUaXDpC/7AI/lZ2DpnUvEO1Sia6WqeigcezDdsRM48zSAbJWTODvXEsn+T2zMaF5+MFZQrE+fLp0gNwVE/hRu/rQ651vwSZ+fAwhvrREpEtndIQuHSB3xQR+6tk9nVkygWmRTZQfnJYZGbQdYZYOkL1iAr8rvDxyKomXdhgPT8ukSwfIXVaBn7qVLsWWNki5rEP8Tt/66+hItPCBuSCrwI/YNWhbT+qscgVb2MENUMaDd16NjpjGSXwga8UEfmMiG5AEim3ir4Mmlmeoa4QuHSB7xQR+/E7b+Dr6oUFTeWgefPwDxzUyEimfLh0gd9ltgBIVDu6EAwmxoi3Uwo8O2o5arEvnyd++qJ9u3Fb738/M9IdHLtGSvRfWrgNQumICP/k69k0sXhYctI3UILqWzshIrEvnGzc+rK/f+HDt8yVp7QmH6q/WvD50DaBkeQV+A6Edns8er0ItTdQ73MIPfeDEbrz63fYd2nf3+fr+X7611vkfuegWbX32xdrlA8Mgq8CP9sOn1MhaNoH//2YGjWPTMiP//zvGXAvnjejg/feodf6SvRfqqd9tr18BYAhkFfgpRZc2aKRPKdLCjg4auyswZqtRk8YC5b805ppfY/P0rv33WKA7Hv2NvnPLL2pfA0hl0R7ztfr1y1ovp7jAn6szA13xxctSjiFEl0feMT6ueaP1K3DEAXvphxse1zlX3F37GkAqRy3bp4zAN7PVkv5B0qikb7r7+f2OPXLp3pFyap8rNbM8sdTdKrDGVRrZ4rD++dFBWwsuj7xj3DUa+BPj06ccqQ+fcOic/cDHcIs0dmZVTpsXN7NRSV+V9E5JmyXdambr3P3eXsfPb+J/eo7+wkdXqzSLdelEp2WOBBdv2zE2rvkj9bt0zExL99mt9vnAMGj7xqvjJG1y94fdfbukSyWtabnMWppaBydynSaWOK59bnADlNERC/Xh7xjzgbVygGHVduAfJOnRCT9vrh5r3M4ulWATP7qnbu21fNzj8/CjG6CEp2XWP/+lcde8wKAtgOm1/RvWK0JeFgtmdoaZrTez9du2bWu5Ov2l7gmKr4NjwS0Wo4O2TXTp0MIH2tR24G+WdPCEn5dLemziAe5+obuvcvdVS5YsqV1QU3fa1l4eeWeXTr3QC6+WmXgDlNER01igiU+XDtC+tgP/VkkrzewwM1sg6XRJ69oscM5ugKJYH3rqDVCi6+G/ND4emocPYHqtztJx9x1m9klJ16ozLfMid9/QZplR9fvwO+r34TfRwq9/viu2AUp0x6sdY655dOkArWp9Hr67Xy3p6rbLCQduQ4O9tf/CCJefdgOU6NIKL42NM2gLtIzfsEx4sBM/vMVidHnkaB/+uDdzHwaAvooJ/Oidtk2VH2llh/4PoksrhO+0jXbpjGs0cOMVgOkVuJZOzVkyOxdPa7AysxQbtI0l/vh4fC2d7WPj+shFt9Q6f8vTL+jYQ2jhA20qJvBTb4DSVX+WUHwMIfLXRbSF/7aVi3Xzw7/SM8+/VOv81y3bRycftbR2+QCmV0zgd8U3EonNw68rvJaO4hugRJz4msW64hOLYxcB0Co6TTOSfnnkTP5MAtCKYgI/GlXhLpXuWjqBG7+iffgpN0ABkL9iAr+rqVUvB62J+wBSrqUDIH/lBH50A5SmlkeuGbvhO20V3VM3NmgLIH/lBH4l3FIOnhdaDz/aiR8QXUsHQP6KC/y6Uq+dFttg8eVbLNYtn0FboGzFBP7OqAp3zQSXR65Zbuodt6IboADIXzGBH5V+eeTgjlPhHbfo0gFKV0zgN7cBSt3zutMy06ylE9+AhUFboHTFBH5X8uWR65bv8TttI+VH19IBkL/iAj8q2Tz81H34cgZtgcIVE/jRO12bEgrcyCyd4PLMTMsEyldM4EftWh45Xew1UXTtD7xglxKA/BUT+KnDyoKd+E3saRvBoC1QvmICv6v20gbBcncNmqYqP9alNU4LHyheMYE/17MqvFpmeC0fWvhA6YoJ/K76fdgNTctMNGgcXcsnugEKgPwVF/gRKRu4Hr3TNnofgFhLByhdMYHfxBaDofKj14nOww/e6csGKED52NN2gtgsmWDgKninbXXu33z/Xi2YN/vP8c2/fl6v2nNB/QoAyF4xgR+5aUlKf8OWFPt/OGrZPlq6z0LdcN/WmmVLbzp0Ue3yAeSvmMDvCi1e1sgsmXqiyxOfeMRi3fy5k+tfAEDxiunDj4ounhYvHwDaVU7gNzAtMrqnbKR81qMH0LZyAj8o3IcfXLyscwkiH0B7ign8JqIy7Tx8WvgA2lVM4KcW3VM3MtgMADNB4FdSx61LNPEBtKqYwN9141PgGqENSDpfU934BQDTKSbwo8JbDEZ33AqulgkA0wkFvpl90Mw2mNm4ma2a9Nw5ZrbJzDaa2amxas6gLtlcpJ7U9wEAKF/0Ttt7JL1f0tcnPmhmR0k6XdLRkl4t6UdmdqS7jwXLm1b9DUgaWh659nr0dOkAaFeohe/u97n7xh5PrZF0qbu/6O4/k7RJ0nGRsgYhdeDSowOgTW314R8k6dEJP2+uHmtNeAOScB9+rPxOC5/EB9Ceabt0zOxHkg7s8dS57n5Vv9N6PNYzCs3sDElnSNIhhxwyXXVa1UQLu/7nDX34ANo1beC7e50lGDdLOnjCz8slPdbn+hdKulCSVq1aVTv1km+AEi2fTcQBtKytLp11kk43s4VmdpiklZJuaamsl4nNgw/Mww/uOAUAbYtOy3yfmW2WdIKkH5jZtZLk7hskfVfSvZJ+KOnMtmfoxDdASb88MvPwAbQpNC3T3a+UdGWf586TdF7k+nXENkAJFBwcNOYPAwBt407bSvrAdeboAGhVMYHfxFo2ofKbuAaJD6BFxQR+E0I9OsHF25ilA6BtBH4ldY9O6vIBlK+4wA8tjxxoYu+807b2WjrOnbYAWlVc4NfV1KBt7AOnmToAQC/FBP6u1nlgWmao/MDJYk9bAO2LLo9cjKaWR/7B3Vu0ZK+Fsz5/y29e0G7zi/n8BZChYgI/9QYoS/baTZL0pWt7rRY9Mye/7oD6FQCAaRQT+F2p7nR968rFuvXck7V9bLz2Ner8ZQAAM1Vc4EdE/0pYsjeBDSBfxXQap77TFgByV0zgAwCmVkzgN7E8MssTAyhZMYHflX7VSwDIUzGB30TjnAY+gJIVE/hdtdeyabgeAJCb4gI/ggY+gJIVE/g7V9Jhi0EA6KmYwG8Cs3QAlKyYwI+vVkkTH0DZign8rtB69M1VAwCyU1zg10UfPoDSFRT41SbikQ1QaOIDKFhBgR9DAx9A6YoJfFrnADC1YgK/KzYPn08NAOUqLvABAL0VE/jxtrnTLQSgaMUEPgBgasUEfhPLItDAB1CyYgK/i8XTAKC3YgK/idY5ffgASlZM4HfV3gCFFj6AwhUX+BHRjdABIGfFBD7LIwPA1EKBb2ZfMrP7zewuM7vSzPab8Nw5ZrbJzDaa2anxqs5MaHlkGvgAChZt4V8v6fXu/geSHpB0jiSZ2VGSTpd0tKTVkv7ZzEaDZU0p3MKngQ+gcKHAd/fr3H1H9eNNkpZX36+RdKm7v+juP5O0SdJxkbJmXKdBFAIAc1CTffgflXRN9f1Bkh6d8Nzm6rFXMLMzzGy9ma3ftm1bg9WZHRc3XgEo27zpDjCzH0k6sMdT57r7VdUx50raIemS7mk9ju/Z+Hb3CyVdKEmrVq2q3UDvzrBx+mYAoKdpA9/dT57qeTNbK+ndkk7yXWm7WdLBEw5bLumxupUcBPdmlmcAgFxFZ+mslvRZSe9x9+cmPLVO0ulmttDMDpO0UtItkbKmr0yrVweAOW/aFv40/knSQknXV63jm9z94+6+wcy+K+ledbp6znT3sWBZM1K3Q4d5+ABKFwp8dz9iiufOk3Re5PoAgOaUc6dt9ZUbrwCgt2ICP4weHQCFKybwG9kAhRY+gIJFB22zc9al/6fd589+FYfHn35B++05v4UaAUAeign8Y5bvpw+8abme275j+oN7WLl0Lx1/+KsarhUA5KOYwN93j/n68gffkLoaAJCtYvrwAQBTI/ABYEgQ+AAwJAh8ABgSBD4ADAkCHwCGBIEPAEOCwAeAIWE5bQloZs9K2pi6HlNYLOnJ1JWYAvWLybl+OddNon5R0fod6u5LpjsotzttN7r7qtSV6MfM1lO/+qhffTnXTaJ+UYOqH106ADAkCHwAGBK5Bf6FqSswDeoXQ/3qy7luEvWLGkj9shq0BQC0J7cWPgCgJUkC38xWm9lGM9tkZmf3eH6hmV1WPX+zma0YYN0ONrOfmNl9ZrbBzM7qccw7zOxpM7uj+u/zg6pfVf4jZnZ3Vfb6Hs+bmX2lev3uMrNjB1i31054Xe4ws2fM7FOTjhno62dmF5nZVjO7Z8Jj+5vZ9Wb2YPV1UZ9z11bHPGhmawdUty+Z2f3Vv92VZrZfn3OnfB+0WL8vmNkvJ/z7ndbn3Cl/z1us32UT6vaImd3R59xBvH498yTZ+8/dB/qfpFFJD0k6XNICSXdKOmrSMZ+Q9LXq+9MlXTbA+i2TdGz1/d6SHuhRv3dI+v6gX7sJ5T8iafEUz58m6RpJJul4STcnqueopMfVmSOc7PWT9HZJx0q6Z8Jjfyvp7Or7syV9scd5+0t6uPq6qPp+0QDqdoqkedX3X+xVt5m8D1qs3xckfWYG//ZT/p63Vb9Jz/+dpM8nfP165kmq91+KFv5xkja5+8Puvl3SpZLWTDpmjaSLq+//XdJJ1sQu5TPg7lvc/fbq+2cl3SfpoEGU3aA1kv7VO26StJ+ZLUtQj5MkPeTuP09Q9k7ufqOkpyY9PPE9drGk9/Y49VRJ17v7U+7+a0nXS1rddt3c/Tp37+7VeZOk5U2WORt9XruZmMnvedhU9asy448lfafpcmdqijxJ8v5LEfgHSXp0ws+b9cpA3XlM9cZ/WtLAN5ytupLeKOnmHk+fYGZ3mtk1Znb0QCsmuaTrzOw2Mzujx/MzeY0H4XT1/2VL+fpJ0lJ33yJ1fiklHdDjmBxex4+q89daL9O9D9r0yarL6aI+3RE5vHZvk/SEuz/Y5/mBvn6T8iTJ+y9F4PdqqU+eKjSTY1plZntJulzSp9z9mUlP365ON8UbJP2jpP8YZN0kvcXdj5X0LklnmtnbJz2fw+u3QNJ7JH2vx9OpX7+ZSvo6mtm5knZIuqTPIdO9D9ryL5JeI+kYSVvU6TaZLPl7UNKHNHXrfmCv3zR50ve0Ho+FXsMUgb9Z0sETfl4u6bF+x5jZPEn7qt6flbWY2Xx1/nEucfcrJj/v7s+4+2+r76+WNN/MFg+qfu7+WPV1q6Qr1fnzeaKZvMZte5ek2939iclPpH79Kk90u7mqr1t7HJPsdawG6N4t6U+96tCdbAbvg1a4+xPuPubu45K+0afcpO/BKjfeL+myfscM6vXrkydJ3n8pAv9WSSvN7LCqFXi6pHWTjlknqTsi/QFJP+73pm9a1e/3LUn3ufsFfY45sDumYGbHqfM6/mpA9dvTzPbufq/OAN89kw5bJ+kj1nG8pKe7fz4OUN/WVcrXb4KJ77G1kq7qccy1kk4xs0VVt8Up1WOtMrPVkj4r6T3u/lyfY2byPmirfhPHg97Xp9yZ/J636WRJ97v75l5PDur1myJP0rz/2hyhnmLk+jR1RqsfknRu9dhfq/MGl6Td1OkK2CTpFkmHD7Bub1Xnz6a7JN1R/XeapI9L+nh1zCclbVBn5sFNkk4cYP0Or8q9s6pD9/WbWD+T9NXq9b1b0qoB//vuoU6A7zvhsWSvnzofPFskvaROq+lj6owJ3SDpwerr/tWxqyR9c8K5H63eh5sk/dmA6rZJnb7b7vuvO2Pt1ZKunup9MKD6/Vv1vrpLneBaNrl+1c+v+D0fRP2qx7/dfb9NODbF69cvT5K8/7jTFgCGBHfaAsCQIPABYEgQ+AAwJAh8ABgSBD4ADAkCHwCGBIEPAEOCwAeAIfH/H6eX5I+9VRkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNX5+PHPmZnsG2RhTSABwp6QsAkiOyqbSK3r1w3XWguItv70q61a2m+LrVWrYq0L4lIV1FYtoiAIIoqyg+xhCWQC2UP2bWbO74+ZhCRMkgEmGTLzvF+vvJKZe+69T+5Mnpw599znKq01QgghvIvB0wEIIYRwP0nuQgjhhSS5CyGEF5LkLoQQXkiSuxBCeCFJ7kII4YUkuQshhBeS5C6EEF5IkrsQQnghk6d2HB0drePj4z21eyGEaJe2bduWp7WOaamdx5J7fHw8W7du9dTuhRCiXVJKHXelnQzLCCGEF5LkLoQQXkiSuxBCeCFJ7kII4YUkuQshhBeS5C6EEF5IkrsQQnghj81zF2fYbJqv9meTnlfGlIGd6R0T6umQhBDtnCR3D9Na85sPd/HvHZkA/O2rQ7z8P0OZMrCzhyMTQrRnMizjYR9uM/PvHZnMn9SH7x6dRP8uYTzwwQ5Onq7wdGhCiHZMkrsHlVdb+PPK/YyMj2TBlL507xDEyzcPxWLT/OXLA54OTwjRjkly96B/b8+ksLyG/ze1HwaDAiC2YzB3XZbAJztPcjS3tNVjKKms4cej+RzLK0Nr3er7E0K0DUnuHvTOpuMkdY9gWM+ODZ6fMyYeP6PinR9cqg903t794Tij//w1N7z6AxOfWc9Nr/1ARkF5q+5TCNE2JLl7SFp2CQezS7hueCxKqQbLOoUFMm1wVz7aZqbKYm2V/f/zmyP89pM9pPbowJtzRvDbGQPYe7KY617ZxIl8SfBCtHeS3D3kiz1ZKAVTB3Vxuvyaod0pqbSw4VCe2/e9Jb2Ap788wIykriy9YyQT+3fi7rG9WP6L0VRarPzi3W1U1rTOPxUhRNuQ5O4hq/dlMbRHRzqFBzpdPqZPNB2C/Vix+6Rb92ux2vjff/9E945BPH1tMkbDmU8NA7qG89z1Kew/VcyLX6e5db9CiLYlyd0Dispr2HuymHGJTd9Mxc9o4IqBnfn6QA4Wq81t+/54u5nDOaU8Pn0goQFnX+YwsX8nZqd047Vvj8nwjBDtmCR3D/jxWD5aw6hekc22G9+3EyWVFnaZT7tlvzab5p/fHCU5NoIrBzV9kdQj0/pjVIq/fXXQLfsVQrQ9uULVA344WkCAyUBKjw7NtrusTzQGBd8cymNYz+b/EbjiuyN5HM0r47kbhpx1Ere+rhFB3DKqB29sPMZvruhHXGTwBe+7sYpqK2v2Z7PpaD7ZRZUE+Bno0ymMif1iSInr0Gx8QoiWSXL3gB+O5jOsZ0cCTMZm20UE+zEkrgMbDuXy0OV9L3i/b286TlSIP9OTurbY9q7LerH0+3Re3XCUP8wefMH7rqW15v3NGTz71SHySqsIDzQR2zGYyhorX+7J4oW1acyN3s48/T4BZSchIhYmPwHJ17stBiF8gST3NlZSWcP+rGIemJzoUvtxiTG8+HUaRRU1RAT5nfd+S6ssrD+Yw+2j41v8pwLQJSKQ38fvY/yOeeid+Sg3JNnyagvz39/Jmv3ZXJIQyd9vTGFUr6i6k7pFFTX89MVrDNv9IgFU2VcqyoD/zrf/LAleCJe1OOaulFqilMpRSu1pYrlSSr2glDqslNqtlBrq/jC9x96TxWgNQ+KaH5KpdUlCJDYN208UXtB+vz2US41Vu16QbPdybsx+hu4qD4U+k2R3Lz+v/ZdVWbj1jc18fSCb380cyAf3jmJMn+gGs3UiAgxcduRZgmoTe62aCvjqifPab5vavRyeGwxPdbB/P5dj1R7X9fS+RbNcOaG6FJjazPJpQKLj617gHxcelvfak1kEwOBuES61T+nRAaNBsS39wpL7mv05RAT5MbzR1bBOWaph1WMYLI2Kl9VUwNqF57xvq03zwAc72HGikJf+Zyh3XZZwZkxdazi1G1Y9Ds8NgrJc5xspOYX15ctg3Z8gczvY3DeDqM6FJqr/zrf/EzzXf4btcd2LYd/t8Z9hG1Ku1BNRSsUDK7TWZw2+KqX+CazXWr/veHwQmKC1PtXcNocPH663bt16PjG3aws+2MEPRwv44bHJrq2wezk5nzxGtC0Pw3kOjVhtmpH/t4bLEqP5+42pzhuV5kDaaji0Co6sg+qSpjc4ZgH0ngQ9RoEpAACbtmGxWbDYLNTYarDYLFi1FZu28Y9v0nh7UzoPTklkdmo3bNqGLj6JPvQFtkNfYis4hjYY0T1GY8veg64swgZoBbXvzlIdSKbqTB8yUNggOBodNxJ6jEZ3TUYb/QHQR9bB9rfs/yRCYmDobdB7YoO6OZqG73mNhqPfwKaXwFJ1ZqkpAEb9iuoeY6mssVJttVFtsX/VWKqxVFVhranAWl3JgF1/xr+mCE3DE8FVxhB2dPoZ2CygLWC1orQFbFaw2X8eUrqRAN3o0wpQrfzYF5ACaND2LStt48xR0Qyo3ocfNWetW4MfB/0HNHpW1fum6Fe9Fz/tZF3lx/6gVMCAdvwT1rX9QGUApUApBpb+gL/NSdyGQPZFTna0PfOllKpbf0DWf/G3nj3VtsYUwsH4W1AG+zoGgwJlRNVuw2Ag/PQBOpu/wGCz1K2njX7kJVxNWaehGAxGlMGAwoAyGlAGIwalwGAkKGsr4XvfQVnr/d5GfyqG3oMtfqy9rcGAAfuvaTAoFAqDAo5tgB8W2zs/wJCqKkJMgXDVC206ZKiU2qa1Ht5iOzck9xXAIq31RsfjtcAjWuuzMrdS6l7svXt69Ogx7Pjx1q2dcjGa/Lf1JESH8vrtLb42Z3o3NfV60H5Brr+Zdi9Hr11IaXEmh4kiY8hdJIyZSZmljLKqUsrzDlB2ajtlOXspKz1FuTJQFhBCWVhnykuzKNMWypWiRiksCixKYUFhVVCDsj82GLCclS6F8H4fZp6if3UNJ4lmpvEVDAoMSmFQCqNBoRQYDQqjOvNz7fK5k/q4NLHBGVeTuztOqDqbs+b0b11r/SrwKth77m7Yd7tSVmXhaF4ZVw3p5toKaxc2TOxwZvy5zxQIjACDEZu2kVOeQ0ZJxpmvjO85kbMLc4SRkshY+7pF78HK987ejwkMHToQbAok2D+MEL8QQvziCClIp4PFgr/WmACTMuLXYzTGjvGYSnMwFWViOp2BqbwAExpTQDimyD6YjP6YzFswWmsw1fY4lQGD1ihtQ4V1xdBjFIYeo1GhnTAoAwqFUgoDBgzHv0Pt+gDK8lAh0aiU/4GEcSgUmYXlbEjLY9ORAkorKhigTjDCdIRhKo2O2n49QO0bsvZ7tTaSozsSpsoJpRKjsp31plWN3o21yzWQFZGCNgWgjYH23rwpEINfAMovCINfEEb/ADr+tBRDdbFj3XobC46CGX8Dg+nMlzKC0QTKBAYDLL8NSrPPfl1CO8MN7zT59gBg2a2eWbfZ9TvBz5cA2v4JRWvAhrbZsGmNzWrFuGIehvKzS2tYgyIpGf97bFZHe5sVre0/a23FZrPR9bvfNZl00pMfBG0Fm3asZ9+/1vZPSr0Pv9nkuoe6XQPYQNvsMdusjscabDb65n/VYN0eNfZPDl3JZ0ZSV6xao7XGatPYtP26Eqs+87Ot3rJg/5YnNVwodyR3MxBX73Es4N5r5r3EgSz7yVRXx9spMtf9WAOcNJnI8DORQQknXk/FbDKR4R+A2WSkqt67zoSim8VCnMVCcmUl3S0Wwmw2QmyaEJuNYFMQIbEjCYkfT3CfywmJiCPQGHj23PLdy2HtQnSRmZM6im973MeNM399dpyFx+HoOjjyNRxdD5VFzn8f/1CYsxK6ptg/8zalx0QY+1vny7rD9YPtUyoP55Ty47ECjuSWsiy3lCUnrmjyD/d4p0shuCPGkGj8wqIIDI8mOKIT/mHREBwJb8+CYidv24g4uH9107HW6tLX+aesiQuhz4zm1x3/hPN1xz8BnZoYRvP0us2u/yR0u6T5dSf+3vm6k/7Q8qfSbS85xuobiYhj6IzHml/3uZVNr3v7P1tYd7DTdVVErFunC7uLO5L7Z8BcpdQHwCVAUUvj7b4qLdten71fl7DmG9ZUwI//pFwpVocE82loCDsCA7DWS4iBykSsKYSeKpDLMBJn1cRV1xBXVU7XihJMxZlNbFzB78xgdGFaZfL1kHw9Cnjp37v5eFsmo/LKiI8OadiuY08YNsf+ZbPCwiicfnirLoNuLiQNFyilSOwcRmLnesfyubgm/vjiiL//4+Y3OKWJZDPZxVk6tQlp7UL7P+VzOT/SHtf15L4nN/FPxZXXylPrekCLY+5KqfeBCUA0kA08CfgBaK1fUfbu3kvYZ9SUA3c4G29v7GI+oVptsXE8v4zeMaF1N9Fwhz+u2Mc7Pxxn38KpDaYA1rFZse18n20b/8wnhnK+Cg2lQkHPmhoml5WTUGMhVhvpOekpoofe2fxVnE30MoiIgwedzmptVnZxJZOeWc8lvaJ44/bhze5bPzcIVe9Tx4Xu22VuOEdx3olOtK0Lea08ta6buPWEamu4WJN7bkkVN7y6iaO5ZYzqFcnSO0YS6Oee8bE73txMVnEVXzwwtuECrTH/9AGf/fgMn6lyMv1MhBgDmdprBrMJY8iPb0KRmUxbFH5XPknnMbe1vLMLTXROvLbhKP+3cj9/viaJm0b2cNpGa82yN/7GrIynCVbVbtu3yy6CPz4hWlNbnlD1Kr/95CdOnq7gnrEJvPbtMZ796hCPTW88pez8HM4tJSXuzDzz8ppyVu98lU/3vstWVYXyg1ER/ZiXfDeTek4myBRkb3jprzmcU8KUZzfwtH8SN7iys+TrqaixUvDZb+mm3HOF6Z2XJbAhLZcnP91LfFQIo3tHNViuteYvqw7yj8MDCB38W2bmvtb2SdYxlCSEr5PkXs/BrBJW7c1mwZREFkzpS35ZNW99n86943oRHRpwQduurLFiLqzgmtRubMnawqf73mV1xnoqsNHTamN+l0u56rLf0SXCeY+4d0wo0aEBbDyczw0jnLdp7CvTeOZXvcBH941mePyFFx4zGhQv3JjKdf/cxJw3N/P4jAHcMCKOAJORjIJy/vj5PlbtzeamkT2Y8bPpoB644H0KIc6PJPd6PtqWgZ9RcfvoeAB+Ob43/96eySc7Mrl7bK/z3/Du5eSteYqHO5Xw6YkIlpxUhNhsTC+vZnavmQyZ8BQqqPkZNEopJvSLYfXeLCxWGyZjyxcXr9qTRXRoAEN7uHBVqos6hviz7N5RLFi2kyc+3cvTXxygQ7A/macrCDAZeGx6f+4Z20uqOgrhYZLcHbTWfL77FOMSY+gYYr/aMbFzGENiI/h058nzT+67l5P2xYPcFRPBaUMEl1RWMq+gjMkxQwm67Q0Ic36bPWcm9e/ER9vMbD9xmpEJzffEK2usrD+Yw9Wp3d16UhggKjSAt+8cybdpeazdn01xpYU+nUL5WWp3unUIcuu+hBDnR5K7w9G8Mk4WVfKrSX0aPH/5wM48s/oQuSVVxISd+9DMkXULuTsmAj+t+STzFL0cFz5gOn5OiR3gssRoTAbFuoM5LSb3VXuzKKu2MuM8r4JriVKKcX1jGNe36btJCSE8R+7E5PD9kXwAxvSObvD8hH6dAPjmUBMFrZpxrOgYd4XaMGh441TOmcQODS5QclV4oB+jekWx8qdTtDTLafnWDOIigxjdK6rZdkII7yTJ3WHTkTy6RQTSM6rhXYcGdg0nOjSADeeY3E8Un+Duz29BK3gjK5t4i6Vhg4jY84pzdmp3jueXN1sCOKOgnO8O53PdsDi3D8kIIdoHSe4OW9MLGZkQedaJQINBcUmvSLYdd73kbkZJBnf+9wZqKgt5vdREgm40+nUBV7VNG9yFID8jy7Y4uUDJYcl3xzAaFNcOO79/IEKI9k+SO5BTXElOSRVJsc5voDG0R0cyT1eQVVTZ4rZOFpu569OfU1lVxGsBfUm8bxM5E/+K2RZtLwcbEXdBF/OEBJi4dlgsn+w4SXbx2fHklVbx/uYTzE6Rk5tC+DJJ7tjvjgQwuFu40+VDHTeybuluSFlFx7nzk9mUVpfxauSl9LvxI/APYUfEFVxW/QK770q3X35/gRfZ3DuuF1atefHrtLOWLfriABar5v6JvS9oH0KI9k2SO2fujjSou/O55oO6ReBvMrC9maGZ7Pw07vrPbIqsFbwaN4uBV79mL+sKHM6x3/iiT6dQt8QbFxnMraN68q8fT/D9kTNlU7/cc4qPtpm5e2wvese4Z19CiPZJpkICe04W0Ss6hNAA54fD32RgQJcw9mcVO12el7WLu7+4lTys/LP/XQwe/VCD5Wk5pXTvEERIE9s/Hw9f2Y8Nh3K59+1tPDZ9AOXVFv666iCpPTqwYIprN98WQngvSe7Yk2/fzs2X4e3fJZyv9mejtW5w0jU//VvuWvtLsg3wSuqvSUm58+ztZ5eS2Nm9PemQABPv3TOKu9/ewmP/+QmAy/pE8/cbU9xW6EwI0X75fHKvsdo4kV/O9MHNX+zTv2sYy7ZmkFtSRafwQAAK93/KPd89ykmTkZcveYqhA649az2rTXMkt5Qxfdw/37xLRCCf/eoy9mcVE2Ay0DsmVC77F0IAktw5nl+Oxabp3Smk2Xb9u9hPtu7PKqFTeCBFW17n3p3PcMLPj8Vjn2ZE7+lO1zMXllNlsZHYqYUbdJwng0ExyNU7OwkhfIbPJ/cjufa7I7V0AjKpYBUb/X9H9/fyKQ4M597IQI4GBPLi+L9xSfzlTa53IMtxMtXNwzJCCNEcSe6O5N6rueS+ezmhqx8i1FBBiVLc1zGQNH9/nk+4lkubSexgn2ZpUDCgi/NplkII0Rp8firkkZwyOocHNDlTBrDf2aemgjKl+GWXTuwP8OfZ7DzGbVve4vb3nSyid0woQW1wt3MhhKjl88k9o6CcnlHNj7dTZMYGzOscw54Af/6ak8eEigqXin/tySxmcBPz54UQorX4fHI3F5YT27GFy/QjYtkREMCWoEAeyS9kSnlF3fPNySutIqu4kkFNXPkqhBCtxaeTe7XFRlZxJbEdg5tvOOm3rAgNJshm4+rSMvtzLhT/qi1rILNZhBBtzaeT+6miCmyaFnvuVSHRrAoJYXKVjUAbFPl3can4V21Zg4HScxdCtDGfni1jLrQPr8S10HP/ducblBgNXDXtRa761I/IEH/eSb6kxe1vO15In06hRAT5uSVeIYRwlU/33M2F5UALPXdLNSvydxCt/BgZN47+XcLr5q43x2bTbE0vYES8+25OLYQQrvLx5F6B0aDoGhHYZJuig5/zTYCJaZ1HYTKYGNgtnNySqhZru6fllFJcaWF4z+bvdSqEEK3Bp5N7RkE5XcIDMRmbPgyrfnoTi1JclXofAMN62nviLdV233q8AIAR8ZLchRBtz6eTu7mwovkhmepyPi8+RG9DMP1jkgD7PVUDTAa2pjef3LccKyAmLIC4SLkbkhCi7bmU3JVSU5VSB5VSh5VSjzpZ3kMptU4ptUMptVsp5byK1kUm83QF3ZtJ7uY9H7A9wI+ZPSbXVVv0NxkYEteBbc303G02zcbDeYzuFSVVGoUQHtFicldKGYHFwDRgIHCTUmpgo2a/BZZrrVOBG4GX3R2ou9lsmpySqmbH2z/f9x4AM4be3+D5EfEd2ZNZRHFljdP19p4sJq+0mon9Y9wXsBBCnANXeu4jgcNa66Na62rgA+DqRm00UDuZOwI46b4QW0d+WTVWm6ZzuPPkritOs6Iyk+F+kXQNa3gl6vi+nbDaNN+l5Tldd93BHJSCcYmS3IUQnuFKcu8OZNR7bHY8V99TwC1KKTOwEpjnluhaUXaxfbZLpzDnyX3vjtdJ9zMxs/dVZy0b2qMDYYEm1h3Mcbruyp9OkRrXgajQAPcFLIQQ58CV5O5s0Fg3enwTsFRrHQtMB95RSp21baXUvUqprUqprbm5uecerRvllNiTe+dw5wl4Rdqn+Gu4POWes5aZjAbG941hzf4cqi22Bsv2nyrmQFYJP0tt/P9PCCHajivJ3QzE1Xscy9nDLncBywG01puAQCC68Ya01q9qrYdrrYfHxHh2yCKrqArA6bBMTUkWX1jzGR/UnfAA53Vhfj40loKyatbuz27w/LItGZgMihnJ3dwftBBCuMiV5L4FSFRKJSil/LGfMP2sUZsTwGQApdQA7Mnds13zFmQXV6IUxISd3XPftOVFCoxGZg64ocn1x/WNoUt4IO/+eLzuubzSKj7YcoJZKd2IDPFvlbiFEMIVLSZ3rbUFmAusAvZjnxWzVym1UCk1y9Hs18A9SqldwPvAHK1146Gbi0pOSSVRIQH4ObmAacXxr4jQMHbgzU2ubzQo7h6bwHeH81mzz957/+uXB6m22PjVxD6tFrcQQrjCpcJhWuuV2E+U1n/uiXo/7wPGuDe01pVdXOV0vL0sL411lHF1eH/8TM33vm8d3ZOPtpl54IMdDO3ZkW/T8rh/Qu8W78cqhBCtzWevUM0urnQ63r5m6wtUGgzMTLqzxW0EmIy8eccIRiREkpZdyq8m9ubXV/RrjXCFEOKc+GzJ3+ziKpJjzz5ZuuLkRmKVgSF9XLvItmtEEEvvGOnu8IQQ4oL4ZM+9xmojv6zqrJ57tnkzPxpqmBmVImUDhBDtmk8m9/zSarQ+e6bMF9sXo5VixrBfeigyIYRwD99M7mX2Oe5R9acras2KvB0kaT/iu4/yUGRCCOEePpncC8vsBb8iQ8703A8d/oKDRs3MLpd6KiwhhHAbn0zuBeXVAESGnLm36Yrdb2DUmqkjF3gqLCGEcBufTO6FZfbk3jHYPixjs1n5vPgQY1QIkZFyAZIQov3zyeSeX1aNUtDBkdy3/PQvcgxwVY/LPRyZEEK4h08m98KyajoE+WE02Kc7rjjwPiE2G+NHXPSVioUQwiU+mdwLyqvp6JgpU1ldxlflGUwxRREU2tnDkQkhhHv4ZnIvrSbSMSSzfsc/KTMoruo9q4W1hBCi/fDJ5F5YXl1XknfFkc/oZLEyfOgvPByVEEK4j08m94Iye3IvKM3mu+o8ZgTFYgwM83RYQgjhNj6X3LXWFDrG3L/c9hIWpZg54CZPhyWEEG7lc8m9pMpCjVUTGezP5yfW0LfGSt+kpm/KIYQQ7ZHPJffaC5gMKoPdtlJmhidCCzflEEKI9sbnknuBI7mfyHkXpTXTku7wcERCCOF+PprcNT+UbGNkjaZLv6s8HZIQQridTyb3sKBDnFT2m3JgMHo6JCGEcDufS+6F5dXERqwlwGZjSqrMbRdCeCefS+65peUUhmcw0WIgNH6sp8MRQohW4XM3yD5RuJ4yo+aqTpeC3CdVCOGlfKvnvns5AUUvEWm1MvrAOti93NMRCSFEq/Cd5L57OcUrHuC7ID+mlpbjV5oN/50vCV4I4ZV8J7mvXchuo41qg2Jyebn9uZoKWLvQs3EJIUQr8J3kXmTGbLKfYkioqWnwvBBCeBuXkrtSaqpS6qBS6rBS6tEm2lyvlNqnlNqrlHrPvWG6QUQsZj8TATYb0VZbg+eFEMLbtJjclVJGYDEwDRgI3KSUGtioTSLwv8AYrfUgYEErxHphJj+B2eRHrMVC3RwZvyCY/IQnoxJCiFbhSs99JHBYa31Ua10NfABc3ajNPcBirXUhgNY6x71hukHy9ZwIDKFbjRWNgog4uOoFSL7e05EJIYTbuTLPvTuQUe+xGbikUZu+AEqp7wAj8JTW+ku3ROgmWmtOKhudaqLYdMsPXNon2tMhCSFEq3EluTu70kc72U4iMAGIBb5VSg3WWp9usCGl7gXuBejRo8c5B3shTlcUUGYAVdOh7ubYQgjhrVwZljEDcfUexwInnbT5VGtdo7U+BhzEnuwb0Fq/qrUerrUeHhMTc74xnxdzzi4Aqqpj6u6fKoQQ3sqV5L4FSFRKJSil/IEbgc8atfkEmAiglIrGPkxz1J2BXqja5F5UHUuHYD8PRyOEEK2rxeSutbYAc4FVwH5gudZ6r1JqoVJqlqPZKiBfKbUPWAc8rLXOb62gz4e58DAApfQmwCRlfoUQ3s2lwmFa65XAykbPPVHvZw085Pi6KJlLMoiyWKkOivd0KEII0ep8piqkuTKfzlZFZWiwp0MRQohW5zPlB8yWUjpYAomSk6lCCB/gE8m9xlpDFlb8LWF0DJbkLoTwfj6R3E8VHsamFDVVHYkMkZkyQgjv5xPJ3Zy1A4CSyi5yAZMQwif4RnLP3w9AfnUPImVYRgjhA3wjuRcdxU9rMmsS5OpUIYRP8I3kXpZFd4uNEsKICpXkLoTwfr6R3KuL6KztST0yJMDD0QghROvz+uSutSZDVxGlwgBkWEYI4RO8PrkXVxZSqiBEReFnVIQH+sxFuUIIH+b1yb12GqRBd6FjsD9KOStPL4QQ3sXrk3tG7k8AVFjjZEhGCOEzvD65mwvSAMipjpeZMkIIn+H9yb3UTKTVytGKzjJTRgjhM7w/uVfkEauN5JRbiZQ7MAkhfITXJ/dMayndTSGUVFqk5y6E8BlendwtNgunsNLFLwqASBlzF0L4CK9O7lkFh7EqRVRgNwC5UYcQwmd4dXI3Z9vnuIcGxANydaoQwnd4d3LP2weAITARkJ67EMJ3eHdyP30Mk9aUGezJXXruQghf4d3JvTyL7lYbOTXBKAUd5EYdQggf4d3Jvfo0sYZA8suq6Rjsj9EgdWWEEL7Bu5O7rYpY/w4UlFXLkIwQwqd4bXIvriykyACxwV3Il+QuhPAxXpvcM7N2AhAbkUBBWbXMlBFC+BSvTe7mnN0AxEYPkGEZIYTPcSm5K6WmKqUOKqUOK6UebabdtUoprZQa7r4Qz4+50F7qt0vnVArLpecuhPAtLSZ3pZQRWAxMAwYCNymlBjppFwbMB350d5Dnw1ySQQerlerAnmgNMWFSNEwI4Ttc6bmPBA5rrY9qrauBD4CrnbT7A/AXoNKN8Z232lK/uWVWQJK7EMK3uJLcuwMZ9R6bHc/VUUqlAnH0TYOxAAAeeklEQVRa6xXNbUgpda9SaqtSamtubu45B3suzJZSYk2h5JVWA5LchRC+xZXk7uzKH123UCkD8Bzw65Y2pLV+VWs9XGs9PCYmxvUoz5HVZuWkshIbGEVuSRUA0aGS3IUQvsOV5G4G4uo9jgVO1nscBgwG1iul0oFRwGeePKmaXXgYi1LEhsWRVyrJXQjhe1xJ7luARKVUglLKH7gR+Kx2oda6SGsdrbWO11rHAz8As7TWW1slYheYs7YDEBuZSG5JFcH+RkICTJ4KRwgh2lyLyV1rbQHmAquA/cByrfVepdRCpdSs1g7wfNSW+u0ek0RuSZWMtwshfI5L3Vmt9UpgZaPnnmii7YQLD+vCmIuOYdSaLl2HkleaJkMyQgif45VXqJrLsuhq1ZiC7SdUYyS5CyF8jHcm9+rTxBrsCT2vVIZlhBC+xzuTu62KWP8Iqi02CstrZFhGCOFzvC65l1acprCu1K99GqT03IUQvsbrkntm9g7AXuo3r8R+dWp0qBQNE0L4Fq9L7uacnwCIjR5Ibqm9zI303IUQvsb7knvhIQBiu6TW67lLchdC+BavS+4ZJRmEWW1ERPcnt1TG3IUQvsnrkru5Io9YDGAwkltSRVigiUA/o6fDEkKINuV1yT3TUkqsMQSA3FK5gEkI4Zu8KrlbbVYylZXYQHs54dySKqJlSEYI4YO8KrnnFh6mRiliw2Ltj6VomBDCR3lVcs845Sj12zERrTVZRZV0CQ/0cFRCCNH2vCq5m/P3AxDXKYniSgsVNVZJ7kIIn+Rdyb3oGAat6dJ1GNnF9guYOkdIchdC+B7vSu5lp+hq0/gFR5JVZE/u0nMXQvgi70ru1aeJVfYTqFnFktyFEL7Lu5K7o9QvQLaj594pXGbLCCF8j9ck9/LKYgocpX7B3nPvGOwnV6cKIXyS1yR3c7ZjGmREAgDZxZV0liEZIYSP8qLkvguA2OgBgL3nLsldCOGrvCa5ZxYeBiC2y1AAsoqq5GSqEMJneU1yN5dkEGqzERHVjxqrjfyyKpnjLoTwWd6T3CtyidUGlNFEbkkVWss0SCGE7zJ5OgB3MVtK6WUMBerNcY+QaZDi4lFTU4PZbKaystLToYh2IDAwkNjYWPz8/M5rfa9I7jZtI1NZGRcUBVB3daqcUBUXE7PZTFhYGPHx8SilPB2OuIhprcnPz8dsNpOQkHBe2/CKYZm8wiNUKUVsaBwAmYUVAMR2CPZkWEI0UFlZSVRUlCR20SKlFFFRURf0Kc+l5K6UmqqUOqiUOqyUetTJ8oeUUvuUUruVUmuVUj3PO6LzYD61DYDYjn0AyDxdQWiAifAgr/hgIryIJHbhqgt9r7SY3JVSRmAxMA0YCNyklBrYqNkOYLjWOhn4CPjLBUV1jsx59lK/sZ2S7I8LK+jeIUj+kIRoJCsrixtvvJHevXszcOBApk+fzqFDh85rW88//zzl5eV1j6dPn87p06fdFaq4QK703EcCh7XWR7XW1cAHwNX1G2it12mta1/lH4BY94bZPHPRUZTWdOs6HLD33Lt3DGrLEIS46Gmt+dnPfsaECRM4cuQI+/bt409/+hPZ2dnntb3GyX3lypV06NDBXeG6lcVi8XQIbc6V5N4dyKj32Ox4ril3AV9cSFDnylx2ii6OUr8AmYXldO8gyV2I+tatW4efnx/33Xdf3XMpKSmMHTsWrTUPP/wwgwcPJikpiWXLlgGwfv16JkyYwLXXXkv//v25+eab0VrzwgsvcPLkSSZOnMjEiRMBiI+PJy8vj/T0dAYMGMA999zDoEGDuOKKK6iosJ8HmzBhAlu3bgUgLy+P+Ph4wH4+4o477iApKYnU1FTWrVsHwNKlS5k7d25dvDNnzmT9+vVYrVbmzJlTF+9zzz131u87Z84cHnroISZOnMgjjzxCWVkZd955JyNGjCA1NZVPP/0UgL179zJy5EhSUlJITk4mLS2N9PR0+vfvz+23305ycjLXXntt3T+ytWvXkpqaSlJSEnfeeSdVVVV1v/+TTz7J0KFDSUpK4sCBAwB88803pKSkkJKSQmpqKiUlJQD89a9/ZcSIESQnJ/Pkk0+64RVuyJVBaWdjG9ppQ6VuAYYD45tYfi9wL0CPHj1cDLFl9Uv9llTWUFxpIVZ67uIi9vv/7mXfyWK3bnNgt3CevGpQk8v37NnDsGHDnC7797//zc6dO9m1axd5eXmMGDGCcePGAbBjxw727t1Lt27dGDNmDN999x3z58/n2WefZd26dURHR5+1vbS0NN5//31ee+01rr/+ej7++GNuueWWJmNbvHgxAD/99BMHDhzgiiuuaHa4aOfOnWRmZrJnzx6AJoeDDh06xJo1azAajTz22GNMmjSJJUuWcPr0aUaOHMmUKVN45ZVXeOCBB7j55puprq7GarWSnZ3NwYMHeeONNxgzZgx33nknL7/8MnPnzmXOnDmsXbuWvn37ctttt/GPf/yDBQsWABAdHc327dt5+eWXeeaZZ3j99dd55plnWLx4MWPGjKG0tJTAwEBWr15NWloamzdvRmvNrFmz2LBhQ90xdwdXeu5mIK7e41jgZONGSqkpwOPALK11lbMNaa1f1VoP11oPj4mJOZ94nQdYr9Rv5ml7D0GGZYRw3caNG7npppswGo107tyZ8ePHs2XLFgBGjhxJbGwsBoOBlJQU0tPTW9xeQkICKSkpAAwbNqzFdTZu3Mitt94KQP/+/enZs2ezyb1Xr14cPXqUefPm8eWXXxIeHu603XXXXYfRaK8Mu3r1ahYtWkRKSgoTJkygsrKSEydOMHr0aP70pz/x9NNPc/z4cYKC7LkjLi6OMWPGAHDLLbewceNGDh48SEJCAn379gXg9ttvZ8OGDXX7u+aaa876nceMGcNDDz3ECy+8wOnTpzGZTKxevZrVq1eTmprK0KFDOXDgAGlpac0eo3PlSs99C5ColEoAMoEbgf+p30AplQr8E5iqtc5xa4QtqKgqJrdeqd/aaZAyLCMuZs31sFvLoEGD+Oijj5wu09rph3EAAgLOXAxoNBpdGr9uvE7tsIzJZMJmswE0mObX1P7rt6+/TseOHdm1axerVq1i8eLFLF++nCVLlpy1fkhISIN9fPzxx/Tr169BmwEDBnDJJZfw+eefc+WVV/L666/Tq1evsyZkKKWaPU71f+/6x+nRRx9lxowZrFy5klGjRrFmzRq01vzv//4vv/jFL5rd3oVoseeutbYAc4FVwH5gudZ6r1JqoVJqlqPZX4FQ4EOl1E6l1GetFnEjJ7MalvqVnrsQzk2aNImqqipee+21uue2bNnCN998w7hx41i2bBlWq5Xc3Fw2bNjAyJEjm91eWFhY3fixq+Lj49m2zT51uf4/mnHjxvGvf/0LsA+lnDhxgn79+hEfH8/OnTux2WxkZGSwefNmwD5eb7PZ+PnPf84f/vAHtm/f3uK+r7zySl588cW6BL1jxw4Ajh49Sq9evZg/fz6zZs1i9+7dAJw4cYJNmzYB8P7773PZZZfRv39/0tPTOXzYXqjwnXfeYfx4p6PQdY4cOUJSUhKPPPIIw4cP58CBA1x55ZUsWbKE0tJSADIzM8nJcW+/2KWJ4FrrlcDKRs89Ue/nKW6N6hyYc+wvRG2p38zCCvxNBqJDpPSAEPUppfjPf/7DggULWLRoEYGBgcTHx/P8888zbtw4Nm3axJAhQ1BK8Ze//IUuXbrUnRR05t5772XatGl07dq17gRoS37zm99w/fXX88477zBp0qS65++//37uu+8+kpKSMJlMLF26lICAAMaMGUNCQgJJSUkMHjyYoUPtVV8zMzO544476nr1f/7zn1vc9+9+9zsWLFhAcnIyWmvi4+NZsWIFy5Yt491338XPz48uXbrwxBNPUFxczIABA3jrrbf4xS9+QWJiIr/85S8JDAzkzTff5LrrrsNisTBixIgGJ6idef7551m3bh1Go5GBAwcybdo0AgIC2L9/P6NHjwYgNDSUd999l06dOrl0HF2hWvqY0VqGDx+ua8+aX4h/rZrHoqz1fDPtfSI7DeZX/9rOvlPFrPvNhAsPUgg32r9/PwMGDPB0GMIF6enpzJw5s+6Erac4e88opbZprYe3tG67Lz9gLskgyGajY1R/++PTFTLeLoTwee0/uVfkEquNKKN9hCmjoJy4SEnuQojzFx8f7/Fe+4Vq/8ndUkqs0X5GvKiihoKyauKjQlpYSwghvFu7Tu7aZsOMldhAe6nf4/llAMRHS3IXQvi2dp3c808fodKgiA2zl7I5ludI7tJzF0L4uHad3M2nHHPcHaV+j+fbaz/0jJI67kII39auk3tG3j4AYjslA5CeV0bXiEAC/YyeDEuIi5bZbObqq68mMTGR3r1788ADD1BdXd3sOqdPn+bll1+ue3zy5Emuvfbac9rvE088wZo1a84r5vpCQ0MveBu+ol0nd3PRUQC6d7UXQzqWXyZDMsJ77F4Ozw2GpzrYv+9efkGb01pzzTXXMHv2bNLS0jh06BClpaU8/vjjza7XOLl369atyTIGTVm4cCFTprT+tY5a6wblCnxZ+07uZafoZNUEBNeeUC0nPlqGZIQX2L0c/jsfijIAbf/+3/kXlOC//vprAgMDueOOOwB7/ZPnnnuOJUuWUF5eztKlS7n66quZOnUq/fr14/e//z1gr41y5MgRUlJSePjhh0lPT2fw4MGAvSTv7Nmzueqqq0hISOCll17i2WefJTU1lVGjRlFQUADYy+9+9NFHbN26ta78bVJSUl39liNHjjB16lSGDRvG2LFj666MPXbsGKNHj2bEiBH87ne/c/p71ZYYvv/++xk6dCgZGRmsXr2a0aNHM3ToUK677rq6y/wfffRRBg4cSHJyMr/5zW/qYrvvvvsYO3Ysffv2ZcWKFUDzZYivueYapk6dSmJiIv/v//0/gCbLEDf1u7W2dn0fOnupX38AisplGqRoR754FLJ+anq5eQtYGxVXramAT+fCtrecr9MlCaYtanKTe/fuPavkb3h4OD169KirlbJ582b27NlDcHAwI0aMYMaMGSxatIg9e/awc+dOgLMqPO7Zs4cdO3ZQWVlJnz59ePrpp9mxYwcPPvggb7/9dl05XIDhw4fXbefhhx9m6tSpgL2UwSuvvEJiYiI//vgj999/P19//TUPPPAAv/zlL7ntttvqygI7c/DgQd58801efvll8vLy+OMf/8iaNWsICQnh6aef5tlnn2Xu3Ln85z//4cCBAyilGpQJTk9P55tvvuHIkSNMnDiRw4cPN1uGeOfOnezYsYOAgAD69evHvHnzyMnJcVqGuKnfrbW17+Ruq2KUv72W9KEcewGjvp3DPBmSEO7ROLG39LwLtNZObz1Z//nLL7+cqCj7J+FrrrmGjRs3Mnv27Ga3O3HiRMLCwggLCyMiIoKrrroKgKSkpLoiXI0tX76c7du3s3r1akpLS/n++++57rrr6pbX3gDju+++4+OPPwbg1ltv5ZFHHnG6vZ49ezJq1CgAfvjhB/bt21dXrre6uprRo0cTHh5OYGAgd999NzNmzGDmzJl1619//fUYDAYSExPp1asXBw4cYOPGjcybNw84uwzx5MmTiYiwlxkfOHAgx48fZ9CgQXVliGfMmMEVV1zR7O/W2tptcq+qKiGnXqnfg1n25J7YWU64iHagmR42YB9jL8o4+/mIOLjj8/Pa5aBBg+oSZa3i4mIyMjLo3bs327Ztc1rmtiX1y/saDIa6xwaDwWl54L179/Lkk0+yYcMGjEYjNpuNDh061PXoG3MlhsalfS+//HLef//9s9pt3ryZtWvX8sEHH/DSSy/V9aDPtbyvszLIzsoQP//8883+bq2p3Y65Z9aV+o0HIC27hBB/o9SVEd5h8hPg1+i97Bdkf/58Nzl5MuXl5bz99tuAfYz417/+NXPmzCE42H6u6quvvqKgoICKigo++eQTxowZc16lfZtSVFTEjTfeyNtvv03tDXvCw8NJSEjgww8/BOzJedeuXYD9RhcffPABQF1J4JaMGjWK7777rm6oqby8vO7kcVFREdOnT+f5559vkHA//PBDbDYbR44c4ejRo/Tr16/JMsRNcVaGuLnfrbW12+RuzrZ/3IuLsldMO5hdQt8uYS79lxfiopd8PVz1gr2njrJ/v+oF+/Pnqbbk74cffkhiYiJ9+/YlMDCQP/3pT3VtLrvsMm699VZSUlL4+c9/zvDhw4mKimLMmDEMHjyYhx9++IJ+rU8++YTjx49zzz331J1YBXvifuONNxgyZAiDBg2qu7/p3//+dxYvXsyIESMoKipyaR8xMTEsXbqUm266ieTkZEaNGsWBAwcoKSlh5syZJCcnM378+Ab3Xe3Xrx/jx49n2rRpvPLKKwQGBnL//fdjtVpJSkrihhtuqCtD3JTMzEwmTJhASkoKc+bMqStD3NTv1trabcnf91bN489Z61k39T2iOg1m6B++4spBXVj082Q3RimE+1zsJX+XLl3K1q1beemllzwdSpuaM2cOM2fOPOe5+23BJ0v+mksyCLRpoqIHkldaTWF5DYlyMlUIIYB2fELVXurXgDIaOZhVCEA/Se5CnLc5c+YwZ84cT4fR5pYuXerpEFpF++25W0qJNdnPkO/OtM8nHdTN+R3QhRDC17TL5N641O/ujCJ6RgXTMcTfw5EJIcTFoV0m94LTx6gwKGJD7aV+d5tPkxzbwcNRCSHExaNdJndzln2WTWzHPuSWVHGyqJIhsREejkoIIS4e7TO5554p9bsl3V6YKLVHR0+GJES78H//938MGjSI5ORkUlJS+PHHH5ttv3TpUk6ePOmWfX/22WcsWtT8lbnp6em89957btmfr2ufyf20vdRvt67D+P5IHiH+RpKl5y5EszZt2sSKFSvYvn07u3fvZs2aNcTFxTW7jjuT+6xZs3j00UebbSPJ3X3aZ3IvzyLGqgkKjuL7I/mMTIjEz9gufxUh2sypU6eIjo6uu8oyOjqabt26AbBt2zbGjx/PsGHDuPLKKzl16lRdid6bb76ZlJQUKioqGmxvwoQJLFiwgEsvvZTBgwezefNmAAoKCpg9e3bd1aG1xcOWLl3K3LlzAfu0y/nz53PppZfSq1evuvrwjz76KN9++y0pKSkNriAV565dznM3VxUSq/wxF5ZzNLeMm0b08HRIQpyTpzc/zYEC99b17h/Zn0dGOq+aCHDFFVewcOFC+vbty5QpU7jhhhsYP348NTU1zJs3j08//ZSYmBiWLVvG448/zpIlS3jppZd45plnGD7c+QWRZWVlfP/992zYsIE777yTPXv28OSTT5Kamsonn3zC119/zW233ea0cNapU6fYuHEjBw4cYNasWVx77bUsWrSIZ555pq6mujh/7TO52yoZERDDyp9OAXDFoM4ejkiIi19oaCjbtm3j22+/Zd26ddxwww0sWrSI4cOHs2fPHi6//HLAXlCsa9euLm3zpptuAmDcuHEUFxdz+vRpNm7cWFd9ctKkSeTn5zutCzN79mwMBgMDBw4kOzvbTb+lqOVScldKTQX+DhiB17XWixotDwDeBoYB+cANWut094ZqV1NVRraj1O9/d51iSGwEPeUGHaKdaa6H3ZqMRiMTJkxgwoQJJCUl8dZbbzFs2DAGDRrEpk2bznl7rpbKdVbQr34RLk/VuPJmLQ5UK6WMwGJgGjAQuEkpNbBRs7uAQq11H+A54Gl3BwrA7uWcfHkYWim6Ht9CwqnPmZXSvVV2JYS3OXjwIGlpaXWPd+7cSc+ePenXrx+5ubl1yb2mpoa9e/cCtFjud9myZQBs3LiRiIgIIiIiGpTKXb9+PdHR0YSHu3b1uDvLC/s6V3ruI4HDWuujAEqpD4CrgX312lwNPOX4+SPgJaWU0u78d+y4p6TZpIFO9Kwo5mn/11HBQ4AEt+1GCG9VWlrKvHnzOH36NCaTiT59+vDqq6/i7+/PRx99xPz58ykqKsJisbBgwQIGDRpUd3/RoKAgNm3aRFBQwxrzHTt25NJLL6W4uJglS5YA8NRTT3HHHXeQnJxMcHAwb73VxG0BnUhOTsZkMjFkyBDmzJnDgw8+6NZj4EtaLPmrlLoWmKq1vtvx+FbgEq313Hpt9jjamB2Pjzja5DW13XMu+eu4M82ysFD+GB3J2hOZdLJa7XWuH9zj+naE8JCLveTvuZowYUKzJ1vFhWvtkr/O7n7R+D+CK21QSt2rlNqqlNqam5vrwq7rKTIDEGO1MrGsnGirtcHzQgghznBlWMYM1L/SIRZofFVDbRuzUsoERAAFjTektX4VeBXsPfdzijQiFooymFRewaTyiobPCyHa3Pr16z0dgmiGKz33LUCiUipBKeUP3Ah81qjNZ8Dtjp+vBb5263g7tMo9JYUQwlu12HPXWluUUnOBVdinQi7RWu9VSi0EtmqtPwPeAN5RSh3G3mO/0e2R1t47cu1C+1BMRKw9sV/APSWFaGtaa7nPr3DJhfaP2+09VIVob44dO0ZYWBhRUVGS4EWztNbk5+dTUlJCQkLD2YCunlBtl1eoCtEexcbGYjabOefJBMInBQYGEht7/ucUJbkL0Ub8/PzO6oUJ0VqklKIQQnghSe5CCOGFJLkLIYQX8thsGaVULnD8PFePBposbeBBEte5kbjO3cUam8R1bi4krp5a65iWGnksuV8IpdRWV6YCtTWJ69xIXOfuYo1N4jo3bRGXDMsIIYQXkuQuhBBeqL0m91c9HUATJK5zI3Gdu4s1Nonr3LR6XO1yzF0IIUTz2mvPXQghRDMu6uSulJqqlDqolDqslHrUyfIApdQyx/IflVLxbRBTnFJqnVJqv1Jqr1LqASdtJiilipRSOx1fbVKXWCmVrpT6ybHPs6qyKbsXHMdrt1JqaBvE1K/ecdiplCpWSi1o1KbNjpdSaolSKsdx97Da5yKVUl8ppdIc3zs2se7tjjZpSqnbnbVxY0x/VUodcLxO/1FKdWhi3WZf81aK7SmlVGa912t6E+s2+/fbCnEtqxdTulJqZxPrtsoxayo3eOz9pbW+KL+wlxc+AvQC/IFdwMBGbe4HXnH8fCOwrA3i6goMdfwcBhxyEtcEYIUHjlk6EN3M8unAF9jvnDUK+NEDr2kW9nm6HjlewDhgKLCn3nN/AR51/Pwo8LST9SKBo47vHR0/d2zFmK4ATI6fn3YWkyuveSvF9hTwGxde62b/ft0dV6PlfwOeaMtj1lRu8NT762LuudfdmFtrXQ3U3pi7vquB2rvvfgRMVq1cS1VrfUprvd3xcwmwH+jemvt0o6uBt7XdD0AHpVTXNtz/ZOCI1vp8L167YFrrDZx9l7D676O3gNlOVr0S+EprXaC1LgS+Aqa2Vkxa69Vaa4vj4Q/Y74DW5po4Xq5w5e+3VeJy5IDrgffdtT8XY2oqN3jk/XUxJ/fuQEa9x2bOTqJ1bRx/CEVAVJtEBziGgVKBH50sHq2U2qWU+kIpNaiNQtLAaqXUNqXUvU6Wu3JMW9ONNP0H54njVauz1voU2P9AgU5O2njy2N2J/ROXMy295q1lrmPIaEkTwwyePF5jgWytdVoTy1v9mDXKDR55f13Myd1tN+ZuDUqpUOBjYIHWurjR4u3Yhx6GAC8Cn7RFTMAYrfVQYBrwK6XUuEbLPXm8/IFZwIdOFnvqeJ0Ljxw7pdTjgAX4VxNNWnrNW8M/gN5ACnAK+xBIYx57rwE30XyvvVWPWQu5ocnVnDx3QcfrYk7u53JjblQzN+Z2N6WUH/YX719a6383Xq61LtZalzp+Xgn4KaWiWzsurfVJx/cc4D/YPxrX58oxbS3TgO1a6+zGCzx1vOrJrh2ecnzPcdKmzY+d46TaTOBm7RiYbcyF19zttNbZWmur1toGvNbEPj3yXnPkgWuAZU21ac1j1kRu8Mj762JO7hfHjbkbcYznvQHs11o/20SbLrVj/0qpkdiPc34rxxWilAqr/Rn7Cbk9jZp9Btym7EYBRbUfF9tAk70pTxyvRuq/j24HPnXSZhVwhVKqo2MY4grHc61CKTUVeASYpbUub6KNK695a8RW/zzNz5rYpyt/v61hCnBAa212trA1j1kzucEz7y93nzF289nn6djPOB8BHnc8txD7Gx4gEPvH/MPAZqBXG8R0GfaPS7uBnY6v6cB9wH2ONnOBvdhnCPwAXNoGcfVy7G+XY9+1x6t+XApY7DiePwHD2+h1DMaerCPqPeeR44X9H8wpoAZ7b+ku7Odp1gJpju+RjrbDgdfrrXun4712GLijlWM6jH0MtvY9VjsrrBuwsrnXvA2O1zuO989u7Imra+PYHI/P+vttzbgczy+tfV/Va9smx6yZ3OCR95dcoSqEEF7oYh6WEUIIcZ4kuQshhBeS5C6EEF5IkrsQQnghSe5CCOGFJLkLIYQXkuQuhBBeSJK7EEJ4of8PYc1Ostdi2oEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotoutput(ucont, uopt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the reasons for the popularity of MPC is how easy it is to change its behaviour using weights in the objective function. Try using this definition instead of the simple one above and see if you can remove the ringing in the controller output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(u, x0=x0):\n",
    "    y = prediction(extend(u))\n",
    "    umag = numpy.abs(u)\n",
    "    constraintpenalty = sum(umag[umag > 2])\n",
    "    movepenalty = sum(numpy.abs(numpy.diff(u)))\n",
    "    strongfinish = numpy.abs(y[-1] - r[-1])\n",
    "    return sum((r - y)**2) + 0*constraintpenalty + 0.1*movepenalty + 0*strongfinish"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}